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Abstract 

This  thesis  studies  the  feasibility  of  applying  electrodynamic  propulsion  to  three  dif¬ 
ferent  orbital  applications:  orbital  plane  change,  Molniya  stationkeeping,  and  rendezvous 
and  docking.  Electrodynamic  propulsion  uses  the  forces  resulting  from  electric  currents 
flowing  through  conductors  as  a  spacecraft  travels  through  the  E^th’s  magnetic  field.  A 
vehicle-independent  expression  for  the  specific  power  required  for  any  maneuver  is  derived 
and  used  to  assess  the  feasibility  of  electrodynamic  propulsion.  Analytical  expressions  for 
the  desired  accelerations  and  combined  current-conductor  vector  for  the  plane  change  and 
Molniya  studies  are  developed  based  on  Lagrange's  planetary  equations.  Solutions  to  the 
forced  Clohessy- Wiltshire  equations  are  developed  to  study  in-plane  rendezvous  and  dock¬ 
ing.  Results  show  electrodynamic  propulsion  can  be  used  to  change  either  inclination  or 
right  ascension  of  the  ascending  node  at  rates  of  approximately  0.4  degrees/day  or  higher 
with  current  spacecraft  specific  power  technology.  Electrodynamic  propulsion  can  also  be 
used  to  negate  the  effects  of  the  Earth's  oblateness  on  a  24  hour  Molniya  orbit  at  90*  incli¬ 
nation.  The  energy  required  for  this  maneuver  is  very  sensitive  to  the  right  ascension  of  the 
ascending  node,  which  determines  the  orientation  of  the  orbital  plane  with  respect  to  the 
magnetic  field.  Specific  power  requirements  of  the  non-optimized  acceleration  and  current- 
conductor  algorithms  are  more  than  required  for  conventional  electric  rocket  propulsion, 
but  can  be  met  by  current  spacecraft  specific  power  technology.  Rjendezvous  and  docking 
are  possible  with  electrodynamic  propulsion,  which  offers  the  advantages  of  allowing  a  soft 
dock  -  relative  velocities  and  accelerations  decay  to  zero  as  the  chase  vehicle  approaches 
the  target  -  and  lack  of  a  thruster  plume  to  impart  momentum  or  contaminate  the  target 
vehicle.  Approaches  along  the  target  velocity  vector  with  no  altitude  change  are  possible 
with  current  spacecraft  specific  power.  Approaches  involving  changes  in  altitude  will  be 
possible  when  modest  improvements  in  spacecraft  power  are  made. 
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ORBITAL  APPLICATIONS  OF 
ELECTRODYNAMIC  PROPULSION 


/.  Background  and  Problem  Statement 

The  concept  of  electrodynamic  propulsion  —  using  the  earth’s  magnetic  field  as  a 
source  of  spacecraft  propulsion  —  has  been  described  by  NASA  (12)  and  studied  in  depth 
by  Lawrence  (8)  and  Spenny  (14).  Lawrence  and  Spenny  created  a  preliminary  design  for  a 
large  spacecraft,  called  the  Precision  Orbital  IVacking  Vehicle  (POTkV),  which  was  capable 
of  precisely  tracking  other  vehicles,  docking,  and  performing  gross  orbital  maneuvers.  The 
POTkV  used  electrodynamic  propulsion  as  its  main  propulsion  source. 

The  fundamental  theory  behind  the  Lawrence  design  is  the  Lorentz  equation,  which 
states  the  force  on  a  current  conductor  in  a  magnetic  field  is  (12:131): 

(1) 

where  £_  is  the  force  exerted  on  the  conductor  by  the  magnetic  field,  i  is  the  current  flowing 
through  the  conductor,  ^  is  a  vector  with  magnitude  equal  to  the  conductor  length  which 
points  in  the  direction  of  positive  current  flow,  and  B.  is  the  geomagnetic  field  vector. 

The  POTkV  conductors  must  make  electrical  contact  with  the  plasma  in  the  earth’s 
ionosphere  to  close  the  current  loop.  Lawrence  placed  plasma-generating  hollow  cathodes 
(12:120-121)  at  the  ends  of  POTkV  conductors  to  make  electrical  contact  with  the  plasma 
because  this  technology  is  the  most  promising  for  high  current  densities  (8:1-2). 

One  of  the  drawbacks  of  the  Lawrence/Spenny  design  is  the  large  power  requirement 
(14:19).  One  way  to  reduce  the  power  requirement  is  to  conduct  the  current  through  some 
medium  with  less  resistance  than  the  ionosphere.  Ladouceur’s  closed  loop  circuit  is  just 
such  a  system  (7).  In  Ladouceur’s  concept,  one  half  of  a  continuous  conductor  loop  is 
shielded  from  the  magnetic  field.  Theoretically,  the  current  inside  the  shielded  half  of  the 
conductor  loop  will  be  unable  to  interact  with  both  the  magnetic  field  of  the  Earth  and 
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the  magnetic  field  of  the  unshielded  half.  The  torque  usually  ^^sociated  with  a  current 
flowing  through  a  loop  in  a  magnetic  field  will  not  exist  (17:204-205);  instead,  there  will 
be  a  non-zero  resultant  force  which  can  be  used  for  thrust. 

The  POTkV  design  can  be  modified  to  incorporate  the  closed  loop  circuit,  but  are 
there  other  feasible  applications  for  electrodynamic  propulsion  and  the  closed  loop  circuit? 
That  is  the  question  this  thesis  attempts  to  answer.  We  will  derive  a  general  equation  for 
the  electrical  power  and  energy  required  to  perform  an  orbital  maneuver.  Power  and  energy 
will  then  be  used  to  evaluate  the  potential  performance  of  the  closed  loop  conductors  for 
several  candidate  applications. 

Several  criteria  can  be  used  to  select  -  or  eliminate  -  potential  applications  for 
electrodynamic  propulsion.  The  first  criterion  is  to  perform  thrusting  at  low  altitude, 
since  the  magnitude  of  ^  drops  off  with  the  inverse  cube  of  the  distance  from  the  center  of 
the  Earth,  as  we  shall  see  later.  The  second  criterion  is  to  use  relatively  small  thrust.  As 
elegant  as  the  Lorentz  equation  is,  the  resulting  forces  for  any  reasonable  current  level  are 
small,  and  there  is  no  value  in  studying  high  thrust  applications.  A  third  criterion  is  for  an 
appUcation  to  exploit  the  unique  feature  of  electrodynamic  propulsion:  during  rendezvous 
and  dnrking  there  is  no  exhaust  plume  to  impart  momentum  to  a  target  vehicle  from  a 
chase  vehicle  (8:3). 

A  recent  survey  of  electric  rocket  propulsion  (thrust  which  is  dependent  upon  the 
expulsion  of  mass)  by  Janson  and  others  identifies  near  and  far  term  applications  (4). 
Several  of  these  applications  meet  the  first  two  criteria,  and  have  been  selected  for  further 
study  in  this  thesis.  These  include  stationkeeping  the  argument  of  perigee  (w)  of  a  90* 
inclination  Molniya  orbit  (4:9)  and  controlling  the  right  ascension  of  the  ascending  node 
(fl)  of  a  sun  synchronous  low  Earth  orbit. 

Other  potential  applications  were  inspired  by  Lawrence  and  Spenny.  The  POTkV 
design  assumed  a  low  inclination  orbit,  and  classical  electric  propulsion  was  used  for  out  of 
plane  thrusts  (8:4-16).  Therefore,  another  potential  application  is  the  control  of  inclination. 
For  this  thesis,  controlling  inclination  of  a  low  Earth  orbit  in  the  vicinity  of  90*  will  be 
studied.  The  methodology  will  be  nearly  identical  to  the  sun  synchronous  fl  control  case 
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discussed  above,  and  we  will  refer  to  the  two  collectively  as  orbital  plane  change.  To  satisfy 
the  third  criterion,  we  will  study  rendezvous  between  two  spacecraft  in  low  altitude,  low 
inclination  orbit. 

For  each  application  we  will  model  the  electrodynamic  accelerations  and  their  re¬ 
sulting  effects  on  the  orbit,  then  compute  the  power  and  energy  required  for  each  mv 
neuver.  Janson  reports  the  current  state  of  the  art  in  spacecraft  specific  power  is  four 
Watts/kilogram  (W/kg)  (4:3),  and  we  will  use  this  as  a  benchmark  in  evaluating  the 
applications  we  study.  We  will  determine  whether  each  application  is  feasible  for  electro¬ 
dynamic  propulsion  based  on  a)  ability  to  perform  the  stated  objective,  b)  the  electrical 
power  required  compared  to  the  four  W/kg  benchmark,  and  c)  the  existence  any  undesir¬ 
able  side  effects. 

We  will  not  attempt  to  optimize  the  implementation  of  electrodynamic  propulsion 
—  our  goal  is  to  show  it  can  be  done,  not  how  it  can  heat  be  done.  Optimization  is  left  as 
a  challenge  for  future  researchers. 

In  this  chapter  we  have  discussed  the  recent  history  of  studies  in  electrodynamic 
propulsion  and  identified  the  potential  new  applications  for  electrodynamic  propulsion  we 
will  study.  Chapter  two  introduces  all  the  necessary  definitions  and  theory  we  will  use  to 
study  each  application.  Chapters  three  through  five  study  the  areas  of  orbital  plane  change, 
polar  Molniya  stationkeeping,  and  rendezvous  respectively,  and  may  be  read  independent 
of  each  other.  The  final  chapter  summarizes  important  results  and  makes  recommendations 
for  further  research.  The  appendices  contain  computer  file  listings,  and  information  about 
potential  conducting  metals  for  the  closed  circuit  conductor. 
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11.  Theory  &  Approach 


In  this  chapter  we  develop  the  theory  necessary  to  study  the  proposed  applications  of 
elect rodynamic  propulsion.  We  will  define  the  various  coordinate  frames  and  the  rotation 
matrices  between  each  frame.  The  unforced  orbital  motion  will  be  described,  and  then  we 
will  introduce  Lagrange’s  planetary  equations  as  a  tool  to  use  in  studying  the  M(dniya  and 
plane  change  applications.  The  Clohessy- Wiltshire  equations  will  be  presented  for  use  in 
the  rendezvous  problem.  We  will  model  the  magnetic  field  and  the  dectrodynamic  forces, 
and  derive  an  expression  for  the  specific  power  required  for  each  maneuver.  Finally,  we 
integrate  all  the  theory  into  an  overall  approach  to  each  study. 

2.1  Coordinate  Systems 

This  section  defines  the  coordinate  systems  used  in  this  thesis. 

2.1.1  Geocentric  Equatorial  Frame.  The  origin  of  the  geocentric  equatorial  frame 
(t)  is  at  the  center  of  the  earth.  The  fundamental  plane  is  the  equator  and  the  tt  unit 
vector  points  in  the  vernal  equinox  direction.  The  tg  unit  vector  points  in  the  direction  of 
the  north  pole.  The  t]  unit  vector  completes  the  right>handed  orthogonal  set  (1:55)  of  this 
Cartesian  frame.  The  i  frame  is  the  inertial  reference  frame  for  this  study. 

2.1.2  Greenwich  Equatorial  Frame.  The  Greenwich  equatorial  frame  {g)  frame 
differs  from  the  t  frame  by  a  single  rotation  about  13.  The  angle  of  rotation  9f  locates  the 
Greenwich  prime  meridian  with  respect  to  ti.  This  angle  is  called  the  sidereal  time  and  is 
calculated  using  the  following  simple  formula  (1:99): 


9f  —  +  u^{t  —  to)  (2) 

where  9^,  is  the  value  of  9g  at  reference  time  to  and  is  the  angular  velocity  of  the 
earth’s  rotation.  Although  9g,  for  a  given  to  is  readily  available  from  a  source  such  as 
The  Astronomical  Almanac  (16),  we  will  use  9f,  =  0  for  all  cases  in  this  study  without 
changing  the  significance  of  the  results. 
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S.l.S  Geomagnetic  Equatorial  Frame.  For  this  study,  we  place  the  origin  of  the 
geomagnetic  equatorial  frame  (b)  at  the  center  of  the  Earth,  and  the  fundamental  plane  is 
the  geomagnetic  equator.  The  ba  unit  vector  points  to  the  geomagnetic  north  pole,  which 
we  locate  at  78.3”  north  and  69”  west  in  the  g  frame.  The  bi  unit  vector  points  to  11.7* 
south  and  69*  west  in  the  g  frame.  The  b^  unit  vector  completes  the  right  hand  coordinate 
system.  This  definition  of  the  b  frame  is  based  on  Tascione’s  description  of  a  dipcde  model 
of  the  Earth’s  magnetic  field  (15:33-34). 

2.1.4  Spherical  Magnetic  Frame.  The  origin  of  the  spherical  magnetic  frame  (rh) 
is  at  the  center  of  the  Earth.  The  m,  unit  vector  is  collinear  with  the  satellite’s  radius 
vector  r.  The  unit  vector  ritA.  is  perpendicular  to  r,  lies  parallel  to  the  geomagnetic 
equator,  and  points  in  the  direction  of  increasing  geomagnetic  longitude.  The  unit  vector 

completes  the  right  handed  coordinate  system,  pointing  in  the  direction  of  increasing 
geomagnetic  latitude. 

2.1.5  Orbital  Frame.  The  origin  of  the  orbital  frame  a  is  at  the  center  of  the 
Earth.  The  fundamental  plane  of  the  orbital  frame  is  the  orbital  plane.  The  unit  vector  hr 
is  fixed  to  r  and  points  radially  outward.  The  unit  vector  09  is  perpendicular  to  the  orbital 
plane  and  collinear  with  the  satellite’s  angular  momentum  vector.  The  a$  unit  vector  lies 
in  the  orbital  plane  perpendicular  to  r  to  complete  the  right  handed  coordinate  system. 
The  d  frame  rotates  with  respect  to  ^he  i  frame  with  angular  velocity  =  va^  where  i/  is 
the  satellite’s  instantaneous  angular  velocity. 

2.1.6  Orbital  Reference  Frame.  The  orbital  reference  i.ame  (c)  is  parallel  to  the 
d  frame  and  the  origin  is  at  the  end  of  r.  If  the  orbit  is  circular,  the  c»  unit  vector  will 
coincide  with  the  satellite’s  velocity  vector.  In  that  special  case,  this  c  frame  is  the  same 
as  Lawrence’s  c  frame  (8:2-1).  The  c  frame  also  rotates  with  respect  to  the  t  frame  with 
angular  velocity  ^  =  vhs. 
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i.2  Rotation  Matrices 


In  this  section  we  will  derive  the  rotation  matrices  we  will  need  to  transform  vectors 
between  the  different  coordinate  systems  defined  in  section  2.1.  We  use  the  notation  72^’, 
where  R}'*  is  a  matrix  which  transforms  a  vector  basis  from  frame  2  to  frame  1. 


2.2.1  Orbital  Frame  to  Geocentric  Equatorial  Frame.  The  orientation  of  a  satel¬ 
lite’s  radius  vector  in  space  at  any  time  is  defined  by  three  parameters: 

•  the  right  ascension  of  the  ascending  node,  H,  the  angle  between  the  vernal  equinox 
and  line  of  nodes. 

•  inclination,  /,  the  angle  between  the  orbital  plane  and  the  equatorial  plane. 

•  argument  of  latitude,  u,  the  sum  of  the  argument  of  perigee  u  and  true  anomaly  i/. 

The  process  of  rotating  from  the  t  frame  to  the  a  frame  is  described  by  rotating  an 
angle  fl  about  the  inertial  13  axis,  then  rotating  an  angle  /  about  the  new  ij  axis,  and 
finally  rotating  an  angle  ti  about  the  newest  is  axis.  The  rotation  matrix  is  formed  by 
the  product  of  three  intermediate  rotation  matrices: 


where 


7^3(n)  = 


cos  n  —  sin  n  0 
sin  n  cos  0  0 

0  0  1 


71,(1)  = 


1  0  0 
0  cos  /  -  sin  / 
0  sin  /  cos  / 


«3(«) 


COS  tt  -  sin  u  0 
sin  u  cos  tt  0 

0  0  1 
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The  resulting  rotation  matrix  is: 


72“  = 


cosHcos  u  -  sinOcos/sinu  -cosHsinu  -  sindcos/cosu  sin  12 sin/ 
sinHcosu  +  cosOcos/sinu  -  sin  12  sin  u  +  cos  12  cos /cos  ti  -cosOsin/ 
sin/sinu  sin/cosu  cos/ 


i.t.2  Geocentric  Equatorial  Frame  to  Greenwich  Equatorial  Frame.  As  men¬ 
tioned  in  section  2.1.2,  the  g  frame  differs  from  the  t  by  a  rotation  about  is  through  the 
angle  $f.  Thus  the  matrix  72**  which  rotates  a  vector  from  the  t  frame  to  the  g  frame  is: 


72*^  = 


coafif  sintff  0 
-sinSg  coaSf  0 
0  0  1 


S.i.S  Greenwich  Equatorial  Frame  to  Geomagnetic  Equatorial  Frame.  Trans¬ 
forming  a  vector  from  the  g  frame  to  the  6  frame  first  requires  a  negative  rotation  about 
gz  through  69*  to  the  intermediate  frame  Then  another  negative  rotation  about  of 
11.7*  completes  the  transformation  to  the  h  frame.  The  result  is: 


72**  = 


cos  69 

-sin  69 

0 

sin  69 

cos  69 

0 

0 

0 

1 

cos  11.7  0  -sin  11.7 
0  1  0 
sin  11.7  0  cosll.7 


cos  11.7 cos 69  —cos  11.7 sin 69  -sin  11.7 
sin  69  cos  69  0 

sin  11.7 cos69  —sin  11.7 sin  69  cosll.7 


2.2.4  Spherical  Geomagnetic  Frame  to  Geomagnetic  Equatorial  Frame.  We  trans¬ 
form  vectors  from  the  m  frame  to  the  h  frame  by  making  a  negative  rotation  about  rh} 
through  the  magnetic  latitude  to  the  intermediate  frame  m*.  Another  negative  rotation 
through  the  magnetic  longitude  about  completes  the  transformation  to  the  b  frame. 
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The  result  is: 


le*"*  = 


COS<Am 

0 

-sin-^m 

0 

1 

0 

sin^m 

0 

casern 

cosA«  -8inA„  0 
sin  cosA„  0 
0  0  1 


coeA*cos0„  -sinA^  -cosA^sin^,^ 
sin  A„  cos  cos  A„  -  sin  A«  sin 

sin^m  0  cos^„ 


We  compute  A^  and  from  the  components  of  r  when  expressed  in  the  b  frame, 
denoted  by  *r,  in  the  following  manner.  By  definition 


T  = 


Then 

*r  =  7l^tn**7V*  *£  =:  rji,  +  f,S,  +  fsi, 

A„,  and  <bm  ue  found  by  taking  the  inverse  trigonometric  functions  of 

tan  Am  =  — 


and 


tan^m  = 


Vri^  +  rj* 


2.S  Orbital  Motion 

In  this  section  we  will  describe  the  Keplerian  motion  of  r  in  the  orbital  plane  (1:72- 
73,  185).  The  magnitude  of  r  is 

1  +  e  cos  V 

where  a  is  the  semi-major  axis  of  the  orbit,  e  is  the  eccentricity,  and  i/  is  the  true  anomaly. 
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Write  r  =  rd,,  treat  a  and  e  as  constants,  and  differentiate  to  find  r  to  be 


(esini/d,  +  (1  +  ecosi/)d«) 


To  find  u  we  first  define  the  mean  motion  (n)  according  to  Kepler's  third  law: 


(4) 


nr  2»  dM 

"  *  r  di 


(5) 


where  n  is  the  gravitational  parameter  of  the  central  body,  and  T  is  the  period  of  the  orbit. 
Integration  of  Equation  (5)  with  respect  to  time  yields  the  mean  anomaly  (M): 


M  =  n(t  -  to)  +  Mo 


(6) 


After  numerically  solving  Kepler’s  transcendental  equation 

M  =:  E  —  eainE 


for  the  eccentric  anomaly  (E),  the  inverse  trigonometric  function  of 


cost/  = 


cobE-'C 
1  -  ecosf 


gives  v. 


(7) 


(8) 


2.4  Lagrange’s  Planetary  Equations 

The  foci  of  the  plane  change  and  Molniya  studies  will  be  how  a  force,  besides  the  two 
body  force,  changes  the  orbital  elements.  For  this,  we  turn  to  the  force  form  of  Lagrange’s 
planetary  equations .  These  equations  describe  the  time  derivatives  of  the  classical  elements 
a,  e,  /,  f),  u;,  and  Af,  due  to  a  perturbing  acceleration,  /,,  defined  in  d  frame  components 
as: 

=  frdr  +  /#d»  +  /ads 
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They  are  (18:37-38): 


da 

di 


2e  sin  1/  *  2o\/l  -  e’ 

— 7,  /'■  ■*  II - /• 

nvl  -  nr 


(9) 


Vl  -  e^  sin  1/  .  y/\  -  e^  /o*(l  -  e*) 

na  '  ^  na*«  V  r 


(10) 


dl  _  rcostt  . 
</<  n«*>/l  —  e’  * 


(11) 


dSl  r  sin  u  . 

dt  na’v'l  —  e*sin/  * 


(12) 


<L; 


\^1  —  COSI/ 


nae 


/f  + 


dM,  _  __1_ 

na  V  o 


VUr? 

nae 


1 

1  -f  e  cos  (/ 


inf! 

nae 


r cot /sin  u  ^ 
no*%/l  —  ^  * 


(13) 


/i  da 


^  o*\f-^di 


(14) 


Unfortunately,  Lagrange’s  planetary  equations  in  this  form  pose  several  problems 
because  of  the  way  we  intend  to  use  them.  Pint,  notice  dM,/dt  has  a  secular  term. 
We  will  be  integrating  over  an  entire  day,  and  this  secular  term  will  eventually  become  a 
nuisance.  Battin  shows  how  to  eliminate  this  problem  from  the  potential  form  of  Lagrange’s 
planetary  equations  (2:487-488),  and  the  same  procedure  will  eliminate  t  in  Equation  (14). 
Recall  M  =  nt  +  Mg  and  differentiate  with  respect  to  time: 


dM 

dt 


dn^  dM, 


(15) 
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Differentiate  n  =  to  find 


dn 

1a 


3n  da 
2a  dt 


(16) 


Substitute  Equations  (14)  and  (16)  into  Equation  (15)  and  simplify: 

l-e» 


dM  1  /'2r  1-e*  \  ,  1  -  \  , 

-T- =  n - ( - cost' I /r - +  — jt)  (^^) 

di  na\a  e  )  not  \  o(l-«’)y 


Thus  we  have  a  new  time-independent  expression  for  the  mean  motion  as  a  function  of  the 
two  body  mean  motion  and  the  perturbing  accelerations. 

Equations  (10),  (13),  and  (17)  pose  a  problem  because  they  are  singular  in  e.  This  is 
undesirable  because  we  intend  to  study  the  orbital  plane  change  of  initially  circular  orbits. 
The  problem  is  solved  for  dM/dt  in  this  manner:  for  our  circular  orbit,  which  is  perturbed 
to  a  near  circular  orbit,  we  can  say  v  ftt  M  and  if  ta  M.  Define  a  new  angle  called  the 
mean  argument  of  latitude  as 

s  =  (j  +  M  (18) 


and  its  derivative  as 


ds  daf  dM 
di~'dt'^~dr 


(19) 


If  we  substitute  Equations  (3),  (13),  and  (17)  into  Equation  (19),  simplify,  and  then  expand 
in  a  Taylor  series  to  first  order  in  e,  we  find  the  ds/dt  term  of  Lagrange’s  planetary  equations 
to  be 


dt 


2  ,  sinu  ,  . 


(20) 


The  last  problem  to  solve,  and  the  most  tedious,  is  the  singuluity  in  e  of  Equa¬ 
tion  (10).  To  do  this  we  introduce  two  of  the  equinoctial  elements  (18:22): 


h  =  esinu; 


k  =  ecosh; 


Note  we  easily  recover  e  from  e  =  +  A*.  If  we  differentiate,  substitute  in  Equations 

(10)  and  (13),  neglect  terms,  and  simplify  the  resulting  expressions  are: 


h 

dt 


esinu; 

de  ,  dL) 

-r-  sinui  +  e  cosw— 
dt  dt 

costt,  _  (r  +  o)8inti hr  ,  hrcot/sinu  , 
fr  Tli  /* - - h 


na 


na* 


na* 


(21) 

(22) 


k 

dt 


ecosu; 

de  du 

—  cosij  —  esintj— 
dt  dt 


sinu,  (r  +  a)cosa  +  ibr  hr  cot /sin  u 
Jr  “T  ~  /#  +  «  /a 

na  na  na* 


(23) 

(24) 


Equations  (9),  (11),  (12),  (17),  (20),  (22),  and  (24)  will  be  used  in  the  plane  change 
study.  The  Molniya  study  will  use  Equations  (9-13)  and  (17). 

2.5  Clohessy- Wiltshire  Equations 

We  will  use  the  Clohessy- Wiltshire  equations  to  study  the  application  of  electro¬ 
dynamic  propulsion  to  rendezvous  and  docking.  The  Clohessy- Wiltshire  equations  are  a 
linearized  set  of  equations  which  describe  the  motion  of  one  body  with  respect  to  another 
body  in  a  nearby  circular  orbit.  We  will  use  the  results  of  Kaplan’s  derivation  (5:105-108). 

Consider  a  body  in  a  circular  orbit  with  the  c  frame  attached  to  it  as  defined  in 
section  2.1.6.  We  will  refer  to  this  orbit  as  the  reference  orbit  or  the  target  orbit,  and  the 
origin  of  the  c  frame  as  the  target.  The  position  of  another  body  orbiting  near  the  target 
p  is  given  by 

p  =  xci  +  ycj  +  zca 
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and  the  relative  translational  equations  of  motion  are: 


X  —  2ny  —  3n*x  = 

(25) 

m 

y  +  2nx  = 

(26) 

z  -1-  n’x  = 

^=/3 

(27) 

where  /r,  /«,  /s  are  the  c  frame  components  of  the  non-gravitational  accelerations  acting 
on  the  vehicle.  When  /^  =  /«  =  /j  =  0,  Equations  (25-27)  can  be  solved  in  closed  form  to 
yield: 


x(0 

»(0 

^(0 


1  .  .  2 
cos  nt  +  —Xo  sin  nt  +  di,  -|-  -tio 
n  n 


y«  -  (3y«  +  6nxa)f  +  +  6xo^  sin  nt  +  ^x,  cos  nt 


Za  cos  nt  -t-  —io  sin  nt 
n 


(28) 

(29) 


(30) 


where  x„  y,,  and  Zg  are  the  components  of  £  at  time  t  =  0. 


2.6  Electrodynamics 

In  this  section  we  will  describe  our  model  of  the  geomagnetic  field,  model  the  forces, 
define  a  new  vector  k  and  derive  an  expression  for  the  specific  power  required  for  any 
maneuver. 


2.6.1  Geomagnetic  Field.  For  this  thesis  we  will  use  Tascione’s  simplified  dipole 
geomagnetic  held  model  (15:33-34).  Tascione  defines  the  magnetic  field  vector  as 


(31) 


where  M  is  the  magnetic  moment,  <!>„  is  the  geomagnetic  latitude  as  defined  by  section 
2.2.4,  and  r  is  the  distance  from  the  center  of  the  Earth.  M  for  the  Earth  is  approximately 
8.05  X  10*‘  kg<meters^/(coulomb*8ec).  To  express  £  in  the  a  frame  we  multiply  by  the 
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appropriate  rotation  matrices: 


•JJ  =  "*£ 


(32) 


S.6.2  Electrodynamic  Forces.  To  determine  the  electrodynauic  force,  we  return 
to  Equation  (1)  and  explore  its  derivation,  which  we  find  in  Orear  (11:350,374).  We  be^n 
with  the  definition  of  force  acting  on  a  charge  q  moving  through  the  magnetic  field  with 
velocity  £: 

£  *  X  £  (33) 


Now  consider  the  differential  force  d£_  due  to  a  differential  charge  dq,  and  replace  v  with 

dlldt: 


Rearrange  terms  and  recognize  the  definition  of  current  is  t  =  dq/dt: 


d£=^dLx&^idlxR 


Finally,  we  assume  the  conductor  £  is  a  straight  line  and  integrate  over  the  appropriate 
limits  to  obtain  Equation  (1) 


The  product  i£  can  be  combined  into  one  vector.  Consider  a  vehicle  with  three 
mutually  perpendicular  conductors  fixed  parallel  to  the  o  frame.  Let  positive  current  tV 
flow  through  conductor  *£  =  {NrLr  0  0}  in  the  plus  dr  direction,  where  Nr  is  the  number 
of  turns  of  a  coiled  conductor  and  L,  is  the  length  of  one  loop  of  the  conductor  exposed 
to  the  £  Held  as  shown  in  Figure  1.  Then 


f  0 


£  =  irLxR={ 


—irNrLrBs 

irNrLrB, 
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Divide  by  the  satellite  mass  m  so  the  acceleration  due  to  the  radial  current/conductor 
becomes 


0 

_  uNtJtt  n 
m 


(34) 


Define  the  parameter 


and  rewrite  Equation  (34)  as 

0 

-lirBa 
KrB, 

Similarly  for  currents  and  t'a  flowing  in  a  positive  direction  through  conductors  ^  =  {0  S$Lt  0} 
and  =  {0  0  N3L3)  respectively,  define 


K,  s 


UN,L, 

m 


and 


«3  = 


iaN^Ls 

m 


Then  the  accelerations  due  to  k«  and  K3  are 


(35) 


(36) 


IS 


Conductor  loop 


L 


Figure  1.  Closed  Loop  Conductor 


Equations  (34-36)  can  be  combined  to  obtain  the  total  acceleration  resulting  from  k,,  K0, 


and  K3: 


K0B3  —  K3B0 


z= 


KaBr  -  KrBs 


-  K0Br 


where  k  is  a  vector  defined  by 


(37) 


K  =  KrOr  -|-  K«d#  -|-  KsOs  (38) 

The  new  nomenclature  k  allows  us  to  refer  to  the  combined  product  of  iNL  in  any 
one  conductor  without  regard  to  the  actual  value  of  either  t,  N,  or  L,  and  retains  the  ‘per 
unit  mass’  definition  we  desire. 

Note  that  although  &  can  be  specified  in  any  direction,  we  cannot  accelerate  in  any 
direction.  Accelerations  resulting  from  kxS_  are  constrained  to  be  perpendicular  to  both 
K  and  ^  due  to  the  nature  of  the  cross  product.  Unless  £  is  perpendicular  to  the  desired 
acceleration  direction,  Equation  (37)  will  generate  unwanted  accelerations  perpendicular 
to  the  desired  acceleration  direction.  We  will  see  the  impact  of  this  phenomenon  in  the 
orbital  plane  change  and  Molniya  stationkeeping  applications. 
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2.6.S  Specific  Power  &  Energy.  Now  we  will  derive  an  expression  for  the  spe¬ 
cific  power  required  to  deliver  the  accelerations  which  result  from  Equation  (37).  The 
total  power  (P)  required  to  drive  a  current  through  the  circuit  which  describes  the  mo¬ 
tor/generator  action  of  electrodynamic  propulsion  is  (8:C-3) 

P  =  iV;  +  t’P  (39) 


where  is  the  voltage  induced  across  the  conductor  as  it  moves  through  the  magnetic 
field  and  R  is  the  total  resistance  of  the  conductor.  Resistance  is  defined  by  (11:337-341) 

R  =  (40) 


where  pj,  is  the  resistivity  in  Ohm-meters,  C  is  the  entire  length  of  conductor  the  current 
flows  through  in  meters,  and  A  is  the  cross  sectional  area  of  the  conductor  in  meters^. 
Note  C  =  2iVI,  assuming  the  length  of  the  portion  of  the  conductor  perpendicular  to  L  is 
neglibly  small  compared  to  L,  because  the  resistance  dissipates  power  as  the  current  flows 
through  the  entire  conductor,  including  the  sUelded  part.  If  we  treat  each  conductor  as  a 
separate  circuit,  the  total  power  is  P  —  P^  +  P«  -f  P3. 

We  return  to  our  definition  of  a  component  of  which  we  now  denote  by  Kj  {j  =  r,  9, 3), 
and  solve  for  ij  to  find 

KiTn 

Substituting  this  result,  Equation  (40),  and  C  =  2NL  into  PR  gives 


=  2«Jm 
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Pm 

NM 


(41) 


Next,  we  assume  a  constant  conductor  cross  sectional  area  A,  and  solve  for  the  con¬ 
ductor  volume  Vj  =  CjAj  =  2NjLjAj  in  terms  of  the  conductor  mass  me^  and  conductor 
mass  density  p<,: 

_  m*  _  me 
“  ejAi  ~  2NjLiAj 
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Replacing  NjLjAj  in  Equation  (41)  gives 

i^Rj  =  (43) 

which  is  the  power  dissipated  by  the  resistance  in  one  conductor. 

The  induced  voltage  term  of  Equation  (39)  is  (pven  by  (8:C-1) 

V;  =  -  (a..,  X  £)  •  (1)  (44) 

Urti  =  V-  vg  h  the  relative  velocity  of  the  conductor  with  respect  to  the  magnetic  field. 
The  velocity  of  the  satellite  in  d  frame  components  is  given  by  Ek)uation  (4).  To  determine 
Vg  we  begin  in  the  g  frame,  where  M.  is  fixed.  The  velocity  of  £  due  to  the  Earth’s  rotation 
then  follows  from 

=  'bio  X  'r 

=  (45) 

where  s  {0  0  .  For  the  rest  of  the  derivation  we  denote  a«i  by 


a«l  =  +  Vr(i,d«  +  Vr,l,aj 


We  proceed  by  expanding  iVi  with  our  new  expression  for  t  and  simplifying  as  follows: 


iVi  = 

=  —i  [^rLr  (VrtliBi  —  VrtlaBg)  +  (VftlaBr  —  VrtIrBa)  +  JVyLa  (VrtIrBa  —  VrtlaBr)] 

=  -  (irNrLr)  -  Vr,laBt)  -  (uNfLa)  (VrtlaBr  -  Vr,lrBa) 

-  {*3^3l>i){VrtlrB0  -  Vrtl$Br) 

=z  -Krtn  (VrtuBa  -  VrtlaBt)  -  K,m  (VrtlaBr  -  VrtIrBa)  -  Katn  (VrtlrB,  -  Vrtl,Br) 

=  -m  («,  (VrtiaBa  ~  VrtlaBa)  +  (Vr,laBr  -  t>r,|,Bs)  +  «3  (VrtIrBt  -  VrtiaBr  )]  (46) 


Finally,  we  substitute  Equations  (43)  and  (46)  into  Equation  (39),  assume 

and  divide  by  the  t  tai  satellite  mass  to  obtain  an  expression  for  specific  power  (P,): 

P$  =  *‘«r  (VrtltBa  -  Vr,l,B$)  —  K#  (VrtIrBt  —  Vr,lrBa) 

-«3  (VrtlrB,  -  Vr,l,Br)  +  A—pnPt  («r*  +  «#*  +  /ts*)  (47) 

m* 

Equation  (47)  provides  insight  to  any  system  we  might  consider.  First,  Equation 
(47)  is  not  configuration  specific,  although  we  must  choose  values  of  the  conductor  mass 
fraction  m/me  to  study.  Second,  we  see  the  icj  terms  will  dominate  P,  if  Kj  is  large,  so  it 
is  important  to  minimize  m/me  «td  the  product  p^Pf  The  limiting  value  of  m/me  is  3, 
which  implies  the  entire  vehicle  is  consists  of  only  the  three  conductors  with  no  bus,  power 
supply,  or  payload.  For  most  of  this  study,  we  will  use  the  value  m/me  =  10,  meaning  30% 
of  the  vehicle  mass  is  in  the  conductors.  To  show  the  affect  of  vuying  m/me,  m/mj  =  6 
(50%  conductor  mass)  and  m/me  =  3  (100%  conductor  mus)  will  also  be  investigated  for 
several  cases.  The  product  p^p^  depends  on  the  conductor  material  we  select.  Appendix  B 
presents  the  results  of  a  survey  of  common  conducting  metals.  Aluminum  has  the  lowest 
PkPc  product,  so  we  use  the  value  {PhPc)m  ~  ^  1^~‘  Ohms’kg/meter’  for  this  study. 
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Although  pjtPc  can  potentially  be  made  quite  small  with  a  superconductor,  this  was  not 
investigated. 

Note  the  last  term  in  Equation  (47)  is  always  positive.  Other  terms  can  be  negative, 
depending  on  the  sign  of  Kj  and  the  sign  of  each  bracketed  quantity.  A  negative  implies 
regeneration,  i.e.  storage  of  energy  rather  than  expenditure. 

For  particular  orbits  some  terms  in  Equation  (47)  are  small  compared  to  others.  For 
example,  consider  a  low  inclination  circular  orbit.  Br  and  B$  are  small  compared  to  ^3, 
and  Vr,i,  »  0  and  a  0.  Therefore  the  second  and  third  bracketed  terms  are  small 
compared  to  the  first.  This  leads  to  the  conclusion  that  the  k,  conductor  is  the  primary 
source  of  any  regeneration  which  may  occur  in  this  special  case. 

In  addition  to  the  instantaneous  power  requirement  given  by  Equation  (47),  we  will 
be  interested  in  the  total  energy  {£)  required  to  perform  a  given  maneuver.  To  find  energy 
we  will  simply  integrate  the  power  over  the  thrust  duration,  or 

^  -  /  P{t)dt 

Jt, 

where  U  ^<1  mark  the  beginning  and  end  of  the  thrust  period.  We  will  perform  this 
integration  numerically  using  the  trapezoidal  rule  (6:979). 

Dividing  the  energy  by  the  thrust  time  gives  the  average  specific  power  7,: 


2. 7  Approach 

In  this  section  we  describe  how  the  theory  we  have  presented  is  applied  to  each  study. 
The  general  approach  is  as  follows. 

1.  Specify  all  constants. 

2.  Specify  initial  values  of  all  orbital  elements. 

3.  Calculate  "r. 

4.  Rotate  r  from  the  d  frame  to  the  b  frame. 
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5.  Calculate  Am  aad 

6.  Rotate  r  from  the  6  frame  to  the  m  frame. 

7.  Calculate 

8.  Rotate  ^  from  the  m  frame  to  the  d  frame. 

9.  Calculate  £.  This  will  be  done  with  algorithms  developed  in  later  chapters  specific 
to  each  case. 

10.  Calculate  and  store  specific  power. 

11.  Calculate  /  =  &  x  JJ. 

12.  Evaluate  Lagrange’s  planetary  equations. 

13.  Numerically  integrate  the  orbital  elements  from  ti  to  tj  with  the  Euler  method 
(6:1063-1064). 

14.  Repeat  steps  3  through  13  until  the  final  time. 

15.  Integrate  specific  power  to  obtain  total  specific  energy. 

16.  Calculate  average  specific  power. 

This  procedure  was  implemented  with  the  Matlab  interactive  software  package  ( 13). 
The  script  and  function  files  are  listed  in  Appendix  A.  Any  case-specific  deviations  or 
simplifications  to  the  procedure  are  described  in  Chapters  3-5. 
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HI.  Orbital  Plane  Change 


In  this  chapter  we  will  investigate  the  applicability  of  electrodynamic  propulsion  to 
orbital  plane  change.  First  we  will  study  inclination  change,  then  we  will  look  at  changing 
the  right  ascension  of  the  ascending  node.  In  both  studies  we  neglect  plane  changes  due 
to  any  other  accelerations,  such  as  Earth  oblateness  or  Sun-Moon  perturbations. 


3.1  Inclination  Change 

Studying  inclination  change  interests  us  for  several  reasons.  First,  Lawrence  re¬ 
stricted  his  studies  to  in-plane  thrusting  and  did  not  study  inclination  change.  Second, 
classical  inclination  change  maneuvers  are  relatively  expensive  in  terms  of  fuel,  so  it  may  be 
advantageous  to  use  electrodynamic  propulsion  rather  than  chemical  or  electrical  thrusters. 


3.1.1  Inclination  Forcing  Function.  To  change  /  we  need  to  find  a  forcing  func¬ 
tion  /a,  to  substitute  into  Elquation  (11)  to  create  a  net  change  in  I.  Recall  Equation  (11) 
is 


dl  __  rcostt 
dt  na^y/i  —  e* 


and  simplifies  to 


£ 

dt 


COSH 

na 


/s/ 


(48) 


for  a  circular  orbit.  Note  a  constant  /a,  is  not  a  satisfactory  forcing  function  because  the 
cos  u  function  in  the  equation  will  result  in  zero  average  /.  Instead,  assume  fa,  is  of  the 
form  A/  cos  u,  where  At  is  a  constant  we  will  derive  shortly,  and  examine  the  average  / 
(/)  over  one  orbit: 


/ 


1  A/  cos*  u  . 

—  /  - Ott 

2w  Ja  na 

A,  /*'  ,  . 

- -  /  cos'  udtt 

2x110  Jo 

At  /«  sin  2u\ I*’ 
2xna  \2  4  /lo 

iL 

2na 


(49) 
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Solving  Equation  (49)  for  Ai  gives 


At  =  2nal 

and  the  forcing  function 

-T 

/a,  =  2mt/coeu 

where  /  is  a  desired  value  we  specify. 

Similarly,  we  can  check  to  ensure  /a,  does  not  change  Q  by  looking  at  the  average  Q 
based  on  Equation  (12): 

dt 


If  we  assume  /  is  a  constant  over  one  orbit,  although  not  strictly  correct,  then  tl  is 

2/ 

(j  s  - — : — ;  /  sin  tt  cos  u  du  =  0 
2ir  sin  I  Jo 

Now  we  must  find  the  vector  k  which  will  result  in  the  desired  force  vector 


rsmu 


sin  u  /-  7  \ 

- : — : ( 2nal cos  u ) 

no  sin  /  '  ' 

2/  . 

— ;  sin  u  cos  u 
sin/ 


L  =  {00fs,} 


To  begin,  find  the  unit  vector  k  as  depicted  in  Figure  2  by  crossing  £  into  and  normal¬ 
izing  the  result: 


k 


BofoOr  —  Brfjkt 


BtOr  -  Bri, 


(50) 

(51) 

(52) 
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Figure  2.  Geometry  of  the  k  Vector 


Thus  we  can  rewrite  £  as 

£  =  « 

Substitute  Equation  (53)  into  the  definition  of  /s  from  Equation  (37),  and  scdve  for  k: 


BtOr  -  Bra, 
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and 


K=  ■  Jh _ 

(Bj  + 


(57) 


Equations  (55-57)  guarantee  achievement  of  the  desired  /s,.  However,  they  will  generally 
also  produce  /,  and  /«  accelerations  as  well.  We  will  see  the  effects  of  these  unwanted 
accelerations  in  the  next  section. 

This  algorithm  haa  another  shortfall.  The  algorithm  makes  no  assumptions  about 
and  there  are  times  in  the  orbit  when  the  denominator  of  Equation  (57)  will  be  near 
0,  resulting  in  unreasonable  values  of  k.  To  avoid  this  problem,  we  will  put  a  limit  on  the 
components  of  k.  Since  the  algorithm  does  not  consider  the  varying  £  held,  there  may  be 
an  inefficient  use  of  power.  To  correct  this  deficiency  crosses  the  line  into  optimization, 
which  is  beyond  the  scope  of  this  thesis. 

S.1.2  Inclination  Change  Results.  Average  inclination  change  rates  of  ±l*/day 
and  ±0.1*/day  at  the  initial  inclination  /,  s  90*  were  studied,  as  well  as  ±0.1*/day  at 
lo  =  45*.  Results  of  each  case  are  summarized  in  Table  1.  m/me  =  10  for  all  cases.  To 
study  the  effects  of  the  rotation  of  the  Si  field*  the  duration  of  each  case  was  one  sidereal 
day  (86 160*).  The  other  initial  orbital  elements  of  each  case  were 

a  =  6  578  km 
h  =  0 
k  =  0 

n  =  0 

tt*  =  0 

In  each  case  we  were  able  to  obtain  the  desired  secular  change  in  /  with  no  net  change 
in  12.  Figures  3  and  4  are  representative  of  the  performance  of  /  and  SI  in  each  case. 

Figures  5  and  6  show  the  case  1  k  and  F,  profiles  respectively,  before  we  imposed 
a  limit  on  k.  After  case  1,  we  imposed  the  limit  that  the  magnitude  of  no  component  of 
K  exceed  200  Amp-m/kg.  Figures  7  and  8  show  F,  and  k  for  case  2.  Note  the  effect  of 
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Table 


Inclination  Change  Study  Results 


Case 

1 

2 

3 

4 

5 

6 

7 

/.  (deg) 

90 

Ei 

90 

90 

45 

45 

Kumit  (A  m/kg) 

n/a 

d]3i 

200 

200 

200 

Target  /  (deg/day 

IDO 

Dl 

BQ 

40 

-.10 

.10 

o 

1 

f  (J/kgxlO«) 

360 

220 

210 

msM 

2.8 

9.6 

42 

26 

24 

.47 

.33 

.46 

1.1 

6736 

287 

269 

79 

100 

229 

202 

Aa  (km) 

1.496 

.349 

.144 

-.061 

.109 

-7.336 

7.404 

Figure  3.  Case  2:  Inclination  Change 
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Figure  4.  Case  2:  Right  Ascension  of  the  Ascending  Node  Change 

limiting  k  in  Figure  7.  Further,  Ks  is  zero,  as  would  be  expected  since  the  numerator  of 
Equation  (55)  gives 


Br 

*  ’ 

0 

B,h, 

B, 

■  X  ' 

0 

>  s  < 

-Brh. 

,  ^3  ^ 

,  /3.  . 

0 

1  / 

Referring  back  to  Figure  3,  note  I  does  not  quite  reach  the  target  value.  This  is  because 
the  limit  on  k  prevented  achievement  of  the  desired  /a,  at  times.  /  can  be  adjusted  to 
make  up  for  the  shortfall. 

In  Table  1  we  see  limiting  k  resulted  in  1/3  less  energy  required  for  case  2  compared 

-T 

to  case  1.  In  cases  3-7  we  changed  the  target  I  to  .l^day.  This  one  order  of  magnitude 
change  resulted  in  two  orders  of  magnitude  change  in  E  and  Ti,  bringing  P,  in  reach  of 
current  spacecraft  specific  power  technology. 

Figure  9  shows  how  the  power  is  divided  between  PR  and  tVj.  There  is  some  negative 
iVi,  but  overall  the  iVi  power  is  negligible  compared  to  the  PR  power,  indicating  the  k? 
terms  dominate  Equation  (47). 

Figure  10  shows  /a,  is  a  cosine  function,  as  desired,  except  for  several  points  near 
revolutions  4  and  12  where  the  limit  on  k  is  invoked,  k  must  be  limited  in  areas  where 
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Figure  5.  Case  1:  ic  Profile  (Amp-m/kg) 
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Figure  7.  Case  2:  k  Profile  (Amp*m/kg) 
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Figure  8.  Case  2:  P,  Profile  (W/kg) 


Figure  9.  Case  2:  i*R  and  tV;  Power 

Br  and  Bt  are  near  zero  compared  to  An  example  of  such  a  point  near  revolution  4 
is  shown  in  Figure  11.  In  these  instances,  S.  i*  nearly  parallel  to  thus  Equation  (57) 
becomes  singular. 

The  unwanted  fr  =  K$Bi  and  /#  s  -‘UrBa  accelerations  rise  and  fall  in  phase  with 
the  relatively  large  B,  oscillations  which  occur  twice  per  day  as  the  geomagnetic  poles 
revolve  through  the  orbital  plane,  as  shown  in  Figure  12.  These  accelerations  manifest 
themselves  in  changes  to  o.  In  cases  1*5,  a  initially  diverges  and  then  returns  to  its  initial 
value,  resulting  in  no  net  change  over  the  24  hour  period  of  interest.  Figure  13  shows  how 
a  varies,  and  is  typical  of  all  the  cases  where  /,  =  90*.  However,  for  /,  =  45*,  a  diverges 
and  never  returns  to  its  initial  value.  Figure  14  shows  o  for  case  6. 

3.2  Right  Ascetuion  of  the  Ascending  Node  Change 

There  are  two  potential  applications  of  AH  maneuvers.  The  first  involves  a  satellite 
at  /  =  90*  whose  line  of  nodes  is  precessed  with  angular  velocity  equal  to  the  Earth's  mean 
angular  velocity  as  it  revolves  around  the  sun,  or  (l»  l*/day. 

The  second  application  is  to  fine  tune  Q,  such  as  might  be  required  to  maintain 
a  specified  angular  separation  between  different  orbital  planes  in  a  constellation.  This 
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Figure  11.  Case  2:  Magnetic  Field  Near  Revolution  4 

would  be  an  expensive  fuel  maneuver  over  the  lifetime  of  a  constellation  of  low  EUirth  orbit 
satellites.  Perhaps  electrodynamic  propuluon  can  perform  small  AH  maneuvers  at  low 
power  using  renewable  electrical  energy. 

3.2. t  Right  Ascention  of  the  Ascending  Node  Forcing  Function.  Development 
of  a  forcing  function  /s^  will  f<dlow  the  same  methodology  as  section  3.1.  Simplifying 
Equation  (12)  for  a  circular  orbit  results  in 

dil  _  sina 
dt  ~  no  sin/''* 

Again  we  find  a  sinusoid  which  eliminates  a  constant  /s  from  contention  as  a  potential 
forcing  function.  Rather,  assume  /a^  is  of  the  form  Bn  sin  u,  where  Bn  is  a  constant  we 
find  in  the  same  manner  as  before: 


5,1  r  “j. 

2w  Jo  na  sin  / 

(58) 

=  ,  /"sin’udu 

2x00  sin //o 

(59) 

/'u  .  M*' 

=  - - r-;  1  r -sinocosu) 

2Tnosin/  \2  /lo 

(60) 
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Figure  12.  Magnetic  Fidd  Components  (Tesla) 
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Figure  13.  Case  2:  Semi-major  Axis  Change 


Figure  14.  Case  6:  Semi-major  Axis  Change 
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Bn 

2na  sin  I 

(61) 

Ba 

=  2naSl  sin  I 

(62) 

fsa 

=  2nansin/sintt 

(63) 

Equation  (48)  is  revisited  to  see  how  affects  /: 


costt  . 

/s 

(64) 

tui 

(2naSlsinIsinit) 

(65) 

no  '  ' 

2nsin/costisintt 

(66) 

Again  we  treat  /  approximately  constant  with  respect  to  u  and  And 


/  = 


20 sin/ 


2ir 


/ 

/  cos  u  sin  u  du  =  0 
Jo 


as  desired. 

The  methodology  used  in  section  3.1  to  And  £  also  applies  here.  DeAne  the  desired 
force  vector  as  ^  =  {0  0  /sq}^.  Then 


k  = 


K  = 

(Bi  +  ff*)*'* 

and  K  =  KK.  As  before  we  will  see  unwanted  accelerations  and  /«,  and  we  must  limit 
the  magnitude  of  k  the  algorithm  might  demand  when  the  denominator  of  k  goes  to  zero. 


3.2.2  Right  Ascension  of  the  Ascending  Node  Change  Results.  Average  0  change 
rates  of  ±l'/day  and  ±0.1*/day  at  initial  inclinations  /,  =  90*  and  I,  =  45*  were  studied. 
Results  of  each  case  are  summarized  in  Table  2.  To  study  the  effects  of  the  rotation  of  the 
£  Aeld,  the  duration  of  each  case  was  86 160*.  The  other  initial  orbital  elements  of  each 


37 


Table  2.  Right  Ascension 


of  the  Ascending  Node  Change  Study  Results 


Case 

1 

2 

3 

4 

5 

6 

lo  (deg) 

90 

45 

45 

45 

45 

90 

Kiimit  (-V  m/kg) 

200 

200 

200 

200 

200 

200 

Target  ft  (deg/day) 

1.0 

1.0 

-1.0 

.10 

-.10 

.10 

£{3 /kg  xlO") 

120 

100 

110 

1.9 

.81 

1.5 

P.  (W/kg) 

14 

12 

11 

.22 

.094 

.17 

(W/kg) 

271 

189 

242 

55 

11 

16 

Aa  (kin) 

-1.59 

5.93 

-9.20 

1.11 

-1.08 

-.531 

case  were 


0  =  6  578  km 
b  =  0 
k  =  0 
ft  =  0 
u,  =  0 

In  each  case  we  were  able  to  obtain  the  desired  secular  change  in  ft  with  no  net 
change  in  I.  Figures  (15  and  (16  are  representative  of  the  performance  of  ft  and  /  in  each 
case. 

Figures  (17  and  (18  show  P,  and  k  for  case  4,  which  is  typical  of  the  other  five  cases. 
In  Table  2  we  again  see  how  the  change  in  ft  from  l.C^/day  to  .I'/day  results  in  two  orders 
of  magnitude  change  in  E  and  bringing  for  the  ft  change  maneuver  below  the  four 
Watt /kg  benchmark. 

In  comparing  cases  1  and  6  in  Table  2  to  cases  2  and  4  in  Table  1  we  see  how  the 
variation  in  the  ^  held  affects  the  total  energy  required  for  plane  change  of  polar  orbits.  In 
Equation  (31)  we  see  B  is  twice  as  strong  at  <l>„  =  90*  than  at  <f>„  =  0.  Since  /3„  =  Bn  sin  u 
while  /a,  =  At  cos  u,  the  /a^  function  is  able  to  take  advantage  of  the  stronger  B  field  near 
the  poles.  The  /a,  function  must  expend  more  energy  to  achieve  the  same  forces  while 
thrusting  over  the  equator. 
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Figure  15.  Case  4:  Right  Ascension  of  the  Ascending  Node  Change 


Figure  16.  Case  4:  Inclination  Change 
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Figure  17.  Case  4:  k  Profile  (Amp-m/kg) 
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Figure  18.  Case  4:  P,  Profile  (W/kg) 
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iVi/m 

i^Rjm 

T. 

(W/kg) 

(W/kg) 

(W/kg) 

•0.0766 

14.1202 

14.0436 

-0.0604 

11.7560 

11.6957 

-0.0487 

9.5775 

9.5288 

•0.0464 

7.5952 

7.5487 

-0.0464 

5.8287 

5.7823 

-0.0522 

4.2549 

4.2026 

2.8389 

2.7809 

-0.0638 

1.6818 

1.6179 
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0.8438 

0.0515 
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0.0371 

2.8319 
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0.0209 

4.2375 

4.2584 

0.0001 

5.8008 

5.8010 

-0.0174 

7.5511 

7.5337 

•0.0325 

9.5160 

9.4835 

•0.0429 

11.6771 

11.6342 

•0.0499 

14.0204 

13.9705 

Figure  19  shows  /sg  is  a  sine  function,  as  desired,  except  for  those  times  when  the 
limit  on  k  is  invoked.  The  unwanted  fr  and  /«  accelerations  due  to  large  B3  components 
again  produce  changes  in  a,  as  shown  in  Figure  20  for  case  4. 


The  il  change  cases  were  chosen  for  investigation  of  the  affect  the  choice  of  m/me 
has  on  7,.  In  addition  to  case  6,  19  other  cases  were  run  at  various  Q  values  and  the  7, 
due  to  the  induced  voltage  and  the  k?  terms  were  tabulated.  Table  3  shows  the  results  of 
these  additional  cases.  All  cases  were  run  at  /  =:  90*  and  m/nte  =  10. 

Figure  21  is  a  plot  of  the  tV^  power,  the  PR  power,  and  the  total  7,  for  all  the 
cases  shown  in  Table  3.  At  H  =  .l*/day,  the  power  (either  expended  or  regenerated)  due 
to  the  induced  voltage  is  about  1/7  the  power  required  to  overcome  the  i^R  losses.  At 
n  =  1.0*/day,  the  iVi  power  is  insignificant  compared  to  the  t^R  power.  The  4  W/kg  limit 
is  exceeded  near  fi  =  .5*/day 
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Figure  19.  Case  4:  Accelerations  (m/s^) 
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Figure  20.  Case  4:  Semi-major  Axis  Change 


Figure  21.  7,  Required  for  ASl  {mfme  =  10) 
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Figure  22.  P,  Requited  for  Ad  (m/me  = 

In  Equation  (47)  we  see  the  power  dissipated  as  heat  can  be  decreased  by  lowering 
m/me  without  changing  the  iVi  power.  Lowering  m/me  implies  more  of  the  total  vehicle 
mass  is  taken  up  by  the  conductor.  In  using  m/nte  »  10  and  assuming  three  orthogonal 
conductors,  30%  of  the  vehicle  is  conductor  mass.  The  conductor  mass  can  be  increased  to 
50%  of  the  vehicle  mass  by  using  m/me  =  0.  If  we  retain  the  three  conductor  configuration, 
the  theoretical  limit  of  m/me  is  3,  meaning  the  conductors  account  for  100%  of  the  vehicle 
mass. 

The  existing  data  can  be  modified  to  study  the  effect  of  differing  m/mc  values.  For 
m/me  =  0)  the  PR  values  of  Table  3  are  multiplied  by  .6  and  added  to  the  unchanged  iVi 
power  to  find  a  new  T,.  Likewise,  multiplying  the  PR  values  by  .3  results  in  a  theoretical 
lower  bound  on  P,.  The  modified  data  is  plotted  in  Figures  22  and  23  for  m/me  =  6  mid 
m/me  =  3  respectively.  Note  the  PR  losses  are  still  significant  compared  to  the  tVj  power. 
The  1.0”/day  goal  is  almost  in  reach  of  the  4  W/kg  power  supply  for  m/me  =  3,  but 
further  m/me  reductions  are  not  possible  without  invalidating  many  of  the  assumptions 
which  went  into  the  derivation  of  Equation  (47).  Although  further  power  reductions  might 
be  found  by  deriving  a  more  general  form  of  Equation  (47),  a  more  promising  way  to  reduce 
7,  is  to  replace  the  aluminum  conductor  with  a  superconductor. 
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Figure  23.  Required  for  Aft  {m/me  =  3) 

3.3  Summary  o/  Orbital  Plane  Change  Reatdts 

In  this  section  we  summarize  the  results  of  the  I  and  SI  change  maneuvers.  We 
have  developed  an  algorithm  which  allows  us  to  use  electrodynamic  propulsion  to  change 
both  I  and  SI  as  desired.  However,  the  maneuvers  can  only  be  achieved  with  reasonable 
power  requirements  when  the  magnitude  of  the  plane  change  is  on  the  order  of  .d’/day. 
Until  spacecraft  power  technology  advances,  this  eliminates  electrodynamic  propulsion 
from  applications  such  as  changing  fl  at  1.0*/<iay  to  produce  a  true  polar  sun  synchronous 
orbit  or  relatively  rapid  A/  maneuvers. 

At  slower  rates  such  as  .10*/day,  A/  maneuvers  are  feasible  for  polar  orbits.  Plane 
changes  at  lower  inclinations  may  be  feasible,  provided  the  change  in  o  can  be  negated. 

Most  of  the  power  required  to  perform  the  maneuvers  is  dissipated  as  heat  due  to 
the  resistance  in  the  conductor.  A  larger  conductor  mass  results  in  less  power,  but  at  the 
expense  of  payload  capability.  Superconductors  should  be  investigated  as  a  means  to  lower 
t^.  v:  power  requirements. 

Additional  study  should  focus  on  constraining  the  change  in  a  and  incorporating  a 
priori  knowledge  of  the  B.  into  the  forcing  functions.  Also,  if  the  only  desire  is  to 
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move  n  from  Hi  to  O3,  less  energy  may  be  required  if  /  and  il  change  simultaneously  so 
the  angular  momentum  vector  follows  a  great  circle  route,  rather  than  simply  precessing 
the  line  of  nodes. 
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IV.  Molniya  Orbit  StatioiUeetping 

Now  we  turn  our  attentioo  to  the  application  of  electrodynamic  propulsion  to  sta¬ 
tionkeeping  the  argument  of  perigee  of  a  polar  Molniya  orbit.  A  Molniya  orbit  is  highly 
eccentric  with  a  12  or  24  hour  period  and  large  inclination.  Thus  if  the  argument  of  perigee 
is  270*,  the  apogee  will  remain  over  the  northern  latitudes  of  the  E^th  for  a  large  portion 
of  the  12  or  24  hour  period.  This  provides  multi-hour  communications  capability  with 
minimum  ground  antenna  tracking  for  those  users  who  do  not  have  visibility  to  equatorial 
geosynchronous  communications  satellites. 

Janson  has  proposed  a  stationkeeping  strategy  for  a  24  hour,  90*  inclination  Molniya 
orbit  using  electric  thrusters  (4:9).  In  this  chapter  we  will  look  at  effects  of  the  Earth’s 
oblateness  on  the  Molniya  orbit,  develop  a  stationkeeping  strategy,  and  compare  our  results 
with  Janson’s. 


4.1  Effects  of  the  Earth's  Obtateness 

The  first  order  approximation  of  the  shape  of  the  Earth  is  a  sphere.  The  second  order 
approximation  is  an  oblate  spheroid,  with  the  polar  diameter  some  21  kilometers  less  than 
the  equatorial  diameter.  Meirovitch  derives  the  components  of  the  disturbing  acceleration 
*/  due  to  oblateness  in  terms  of  the  Earth’s  principal  axis  moments  of  inertia  (9:457-459). 
They  are 

ft  — - - -  (1-3  sin*  ti  sin  /) 

,  —3G  {Ce  —  Ab)  .  n  -Jr 

f$  = - - ^sin2usin*/  (67) 

g  —3G{Ce  —  Ab)  .  . 

73  = - - sin  u  sin  2/ 

where  G  is  the  gravitational  constant,  Cb  is  the  Earth’s  moment  of  inertia  about  the  polar 
axis,  and  Ab  is  the  Earth’s  moment  of  inertia  about  the  orthogonal  axes.  If  we  substitute 
Equations  (67)  into  Equation  (13),  which  we  repeat  here  for  reference. 


Vl  -  e*  cos  1/  ,  Vl  -  e* 


(1  \  ,  rcot/sinu  , 

^  -naWl-e^'’ 
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Figure  24.  Effect  of  Earth  Oblateaess  on  24  Hour  Polar  Molniya  Orbit 
we  can  see  how  the  Earth’s  oblateness  affects  ui.  Our  reference  orbit  is  defined  by 

T  =  86160* 

e  =  .82501 
ut  =  270* 

Figure  24  shows  Au;,uaic  <lue  to  the  Earth’s  oblateness  is  approximately  —.07*  after 
one  revolution.  Left  unchecked,  u  will  continue  regressing  until  the  apogee  is  no  longer 
over  the  North  pole,  making  the  orbit  unusable  for  polar  communications.  The  goal  is  to 
drive  to  0. 

4.2  Desired  Forcing  Function 

Finding  a  suitable  forcing  function  to  produce  the  desired  u  is  not  as  straightforward 
as  it  was  for  the  plane  change  case.  Equation  (13)  contains  all  three  components  of  £, 
so  we  have  one  equation  and  three  unknowns.  Even  if  we  choose  to  neglect  /),  which  is 
reasonable  since  the  cot  I  term  negates  its  effect  on  w  for  a  90*  inclination  orbit,  we  still 
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must  find  some  combination  of  fr  and  ft  to  produce  the  desired  u.  To  solve  this  problem 
we  call  upon  the  method  of  Lagrange  multipliers  (9:54-55),  in  which  we  minimize  one 
function  subject  to  the  constraint  imposed  by  another  function.  From  our  experience  with 
the  plane  change  case  we  know  the  k?  terms  dominate  Equation  (47),  so  a  function  worth 


numnuzing  is 


W  (fi)  s  -I-  icj  -1- 


The  constraint  we  impose  on  H  comes  from  Equation  (13).  We  require  Equation  ( 13) 
to  produce  during  the  thrust  time  At,  and  we  define 


We  write  the  constraint  function  g  as 


where 


gin)  =  -*1  cost//,  -f-  kik, sin  1/ ft  -  —  -  w  =  0 

na*vl  -  e* 


kj  =  1  -J- 


l  +  ecoei/ 


We  simplify  g  by  neglecting  the  third  term,  since  cot  /  =  cot  90  =  0.  Thus 

g  (fi)  =  -*i  cos  I//,  -I-  kik2  sin  i//*  -  u>  =  0 

Last,  we  substitute  in  the  definitions  of  /,  and  /#  in  terms  of  k  from  Equation  (37): 

/,  =  litBs  —  KiBt 
f  -isBr  -  lirBs 

The  result  is 


g  =  -ki  cosi/^KtBs  -  KtBt)  +  kikj  sin  i/ («3.flr  -  «,^3)  -  w  =  0 
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Augment  Equation  (68)  by  multiplying  ^9  times  the  undetennined  Lagrange  multiplier 
A  so  that 


T?  =  W  +  Aj 

7?  =  K*  +  4-  «3  +  A  (-*1  cos  u  {KtBa  -  «3^#)  +  sin  -u)  (71) 


Take  the  partial  derivatives  of  with  respect  to  k,,  k#,  and  k^,  set  each  equal  to 
zero,  and  solve  for  A: 


^  -  A*,*j8ini/53  =  0  ^  =  I7kS?7fl; 

^  =  2«3  -  A*i  cosi/Bs  =  0  ^  A,  =  T.^Bi 

^=:2K3  +  X(kiCO8uB,  +  kikiSiiiuBt)=:0  ^ 


Since  A  =  Ai  =  Aj  =  A3,  we  can  solve  for  k,  and  K3  in  terms  of  «#: 


Kr  =  katan  v  k# 


(72) 


«3  = - (cosvB*  +  kjsinuB,) 

CQBVBi 

Finally,  substitute  k,  and  K3  into  Equation  (70)  and  find  ic#  is 


(73) 


«# 


_ -HfBz  cos  1/ _ 

B3  (*i*3  v  +  ki  cos*  u)  +  ki  (cos  1/B0  +  k^  sin  i/Br)^ 


(74) 


Substituting  Equation  (74)  back  into  Equations  (72)  and  (73)  yields  the  complete  solution 
for  K  required  to  obtained  a  desired  u>. 

As  in  Chapter  3,  we  have  a  completely  analytical  algorithm  to  determine  a  required 

to  meet  the  stated  objective.  However,  it  has  deficiencies  also.  Notice  Equations  (72)  and 

(73)  will  call  for  infinite  k,  and  K3  at  1/  =  ±90*.  Like  the  plane  change  algorithm,  this 

algorithm  makes  no  assumptions  about  the  £  field.  Examine  Equation  (74)  at  1/  =  0  and 

simplify  to  _ 

—1JB3 

k^(Bi  +  Bj) 

In  Figure  25  we  observe  B*  and  B3  both  pass  through  zero  at  i/  =  0,  making  singular 
at  that  point. 
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Figure  25.  Geomagnetic  Field  Strength  for  24  Hour  Polar  Molniya  Orbit,  fl,  =  141" 

We  compensate  for  the  first  problem  by  somewhat  arbitrarily  choosing  to  restrict  our 
thrusting  to  the  region  where  coev  >  .1,  or  -84.3*  <  v  <  84.3*.  This  has  the  added  benefit 
of  not  thrusting  at  higher  altitudes  where  the  magnitude  of  M.  drops  off  with  the  inverse 
cube  of  r,  as  shown  in  Equation  (31).  Another  reason  to  not  thrust  over  the  whole  orbit 
is  the  dipole  model  becomes  suspect  above  2000  km  (10:2>21).  Additionally,  we  will  limit 
the  magnitude  of  any  component  of  k  to  500  Amp'm/kg  throughout  the  thrust  period. 
Thrusting  near  i/  =  0  will  not  be  restricted,  but  the  magnitude  will  be  subject  to  the 
above  limit. 

Earlier  we  defined 

fj  s  ' 

At 

where  At  is  a  thrust  time  which  can  now  be  specified.  To  determine  At,  we  find  E  and  M 
at  cosi/  s  .1  from  Equations  (8)  and  (7),  set  M,  =  —  A/  at  t,  and  solve  Equation  (6)  for 
At: 


cosE 

E 


e  + cost/ 
1  +  e  cos  1/ 
.5662  rad 
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M  =  £  -  e  sin  £ 
=  .4942  rad 


M 

Al 


nAt  +  M, 
2M 
n 

3227* 


We  make  two  more  points  before  analyzing  the  results.  Even  though  our  intent  is  to 
minimize  H,  the  method  of  Lagrange  multipliers  only  guarantees  we  will  find  an  extremum, 
not  necessarily  a  minimum.  To  ensure  H  is  minimized,  we  must  look  at  a  point  near  the 
we  solve  for  and  verify  <  H  (k.^j  +  S&sci)'  Also,  we  only  chose  to  use  a  Lagrange 

multiplier  solution  to  aid  in  solving  for  the  three  unknown  components  of  £  with  respect 
to  one  equation  for  u/.  The  fact  that  we  somewhat  optimized  P,,  which  is  supposed  to  be 
beyond  the  scope  of  the  thesis,  is  a  side  benefit. 

^.3  Results 

The  energy  required  to  stationkeep  u;  depends  on  the  initial  right  ascension  of  the 
ascending  node  (H,).  We  see  this  in  Figure  26,  where  we  plot  total  energy  versus  initial 
n.  The  minimum  energy  is  2.4  x  10’  J/kg  at  12,  =  214*.  There  are  several  local  minima 
with  energies  from  2.6  x  10*  to  2.9  x  10’  J/kg.  Local  maxima  of  5.3  x  10’  J/kg  exist  at 
Qo  =  S9*  and  Q,  =  270*.  The  wide  variation  in  energy  is  caused  by  the  11.5*  tilt  of  the 
magnetic  field.  Figure  27  shows  the  S.  fi«^d  for  fl,  =  89*.  Note  the  relative  flatness  of 
Bi  in  Figure  27  compared  to  B3  in  Figure  25.  Since  £3  couples  into  both  f,  and  /«  in 
Equation  (37),  we  can  expect  large  values  of  and  Kr  are  required  to  generate  the  desired 
fr  and  /#.  Hence  the  larger  power  requirement  where  12,  results  in  a  near  zero  £3. 

There  is  an  unwanted  side  effect  on  a.  Figure  28  shows  how  the  forcing  function 
affects  a  as  a  function  of  12,.  a  changes  by  as  much  as  ±200  km,  which  is  undesirable. 
There  are  several  places  where  Aa  is  0.  Unfortunately,  12,  =  214*  is  not  one  of  them. 
However,  An  does  cross  zero  at  12,  =  141’,  which  is  near  a  local  energy  minimum  at 
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Figure  26.  Specific  Energy  Required  for  Molniya  Stationkeeping 
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Figure  27.  Geomagnetic  Field  Strength  for  24  Hour  Polar  Molniya  Orbit,  Qo  =  89' 


Figure  28.  Change  in  Semi-major  Axis  Due  to  Molniya  Stationkeeping 

Do  =  ISS**.  We  choose  to  tradeoff  an  increase  in  energy  for  no  change  in  a,  and  select  the 
point  n*  =  141*  for  closer  study. 

Figure  29  shows  how  u;  varies  over  one  orbit  at  D,  s  141*.  The  algorithm  produces 
Au  =  -.0013*,  and  we  could  make  Au;  smaller  by  fine  tuning  ui'.  Figures  30  and  31 
show  changes  in  I  and  Q,  We  can  see  what  causes  them  by  looking  at  Figure  32,  which 
shows  a  large  fs,  fz  doesn’t  contribute  to  u;  because  I  =  90*,  but  it  does  change  /  and 
n.  Although  unexpected,  these  changes  are  acceptable  within  the  scope  of  this  analysis. 
Figures  33  and  34  show  that  by  choosing  O,  =  141*  we  are  able  to  keep  Ao  and  Ae  near 
0,  as  desired.  Finer  selection  of  Do  should  drive  Aa  and  Ae  closer  to  0. 

The  total  specific  energy  required  is  2.94  x  10'  J/kg.  Figures  35  and  36  show  the  k 
and  Pt  profiles  respectively.  Peak  power  is  approximately  1000  W/kg,  and  occurs  due  to 
the  large  K3  at  the  beginning  of  the  thrust  period,  ns  would  be  even  larger  if  not  for  a  limit 
of  =  500  Ampere-meters/kg.  If  we  average  the  energy  over  At  the  average  power 

7,  is  91  W/kg,  while  averaging  the  energy  over  the  entire  orbital  period  T  gives  7,  =  3.4 
W/kg. 

Although  restricting  the  beginning  of  the  thrust  period  to  t/  =  -84.3*  solved  the  cos  </ 
singularity  in  Ks,  it  neglects  the  fact  that  B3  also  goes  through  zero  at  =  84.0*,  which  is 
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Figure  29.  Change  in  Argument  of  Perigee  for  Molniya  Stationkeeping 


Figure  30.  Change  in  I  Due  to  Molniya  Stationkeeping 
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Figure  31.  Change  in  H  Due  to  Molniya  Stationkeeping 


Figure  32.  Electrodynamic  Accelerations  for  Molniya  Stationkeeping 
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Figure  33.  Change  in  a  for  Molniya  Stationkeeping 


Figure  34.  Change  in  e  for  Molniya  Stationkeeping 
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Figure  35.  k  Profile  for  Molniya  Stationkeeping 


Figure  36.  Power  Profile  for  Molniya  Stationkeeping 
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also  shown  in  Figure  25.  The  peak  near  1/  =  0  due  to  the  weak  B$  and  components  is 
expected.  S3  does  not  pass  through  zero  again  until  v  =  98*,  hence  the  asymmetry  in  the 
power  profile.  The  k  and  power  levels  are  relatively  small  over  the  range  20*  <  |i/|  <  60* 
where  B$  and  S,  are  relatively  large.  This  highlights  the  necessity  to  develop  algorithms 
which  incorporate  some  prior  knowledge  of  the  B.  field.  There  is  also  some  regenerative 
power  in  the  each  conductor  at  different  times  during  the  thrust  period.  However,  it  is 
negligible  compared  to  the  total  energy  required. 

The  check  on  H  (k,^)  <  H  verified  H  had  in  fact  been  minimized  at 

every  point.  Although  power  is  minimized  at  every  point,  this  does  not  imply  total  energy  is 
minimized.  For  instance,  we  can  look  closely  at  Equation  (13)  and  see  any  energy  expended 
near  1/  =  0  is  wasted.  First,  /s  has  no  affect  on  u  due  to  the  cot  /  term.  Likewise,  the  sin  u 
term  negates  the  affect  of  /«.  An  fr  would  be  useful,  but  we  have  already  seen  B$  and  B3 
are  zero  at  1/  =  0. 

Table  4  shows  how  these  results  compare  with  Janson's  stationkeeping  strategy  for 
several  conventional  electric  thrusters  and  a  monopropellant  system  given  a  four  year 
mission  duration  (4:9).  Janson's  approach  to  stationkeeping  is  to  use  a  constant  +/r  for 
six  hours  per  day  centered  on  apogee.  Oectrodynamic  propulsion  offers  no  advantage 
to  electrical  thrusters  in  terms  of  energy  or  power  consumed.  However,  if  the  energy 
can  be  acquired  from  the  sun  and  stored  on  board  during  the  non-thrusting  phase  then 
no  consumable  fuel  is  required.  If  we  treat  the  mass  of  the  conductor  as  the  analog 
of  propellant  mass  for  comparison  purposes  only,  electrodynamic  propulsion  (EDP)  can 
perform  the  mission  for  an  unlimited  time  with  only  30%  ‘propellant'  mass,  compared  to 
the  38%  required  for  a  stationary  plasma  thruster  (SPT). 

There  au'e  several  other  advantages  to  using  electrodynamic  propulsion  in  this  sceanrio 
compared  to  Janson's.  Since  all  the  thrusting  is  done  near  perigee  when  the  payload  would 
not  be  required  for  communications,  strict  attitude  requirements  which  might  not  be  met 
because  of  unwanted  torques  could  be  relaxed.  Electrical  power  normally  allocated  to 
the  payload  could  be  diverted  to  the  conductors.  Any  concerns  about  the  possible  adverse 
affects  of  the  presence  of  a  electrical  thruster  plume  in  the  line  of  sight  between  the  satellite 
and  Earth  would  also  be  eliminated. 
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Table  4.  Specific  Power  Requirements  for  a  24  I: 

our  Polar  Molniya 

Satellite 

Thruster 

Monopropellant 

Arcjet 

SPT 

Ion 

EDP 

Energy  (J/kg) 

60,480 

82,080 

116,640 

294,000 

P,  peak  (W/kg) 

- 

2.8 

3.8 

5.4 

997 

laaEsaiHi 

- 

2.8 

3.8 

5.4 

91 

7.  {£/T) 

- 

.70 

.95 

1.35 

3.4 

Duty  cycle  {At IT) 

- 

.25 

.25 

.25 

Propellant  fraction 

98^ 

38% 

23% 

- 

Even  with  the  deficiencies  in  the  forcing  function  algorithm  which  cause  excessive 
energy  and  power  requirements,  the  7,  is  feasible  today.  Incorporating  knowledge  of  sin¬ 
gularities  and  peaks  in  the  field  into  the  power,  energy,  and  u;  functions  and  approaching 
the  problem  with  the  goal  of  minimizing  total  energy  should  give  electrodynamic  propulsion 
a  clear  P,  advantage  over  electric  rocket  propulsion. 
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V.  Rendezvous  &  Docking 


In  this  chapter  we  will  study  the  application  of  electrodynamic  propulsion  to  ren¬ 
dezvous  and  docking.  A  unique  feature  of  electrodynamic  propulsion  is  there  is  no  exhaust 
plume  to  impart  momentum  on  the  target  as  the  chase  vehicle  approaches.  Additionally, 
if  we  presume  the  electrical  power  comes  from  some  renewable  source,  we  should  be  able 
to  thrust  continuously  and  affect  a  soft  dock  such  that  the  relative  velocity  and  accder- 
ations  decay  to  zero  as  the  chase  vehicle  approaches  the  target.  This  is  in  contrast  to 
a  conventional  two  bum  maneuver,  in  which  the  chase  vehicle  arrives  at  the  target  with 
non-zero  velocity  and  thraster  firings  are  required  to  change  the  relative  velocities  to  zero 
near-instantaneously. 

For  simplicity,  we  take  the  11.7*  tilt  out  of  the  magnetic  field  model  and  study  a 
zero  inclination  reference  orbit.  In  this  scenario  the  M.  field  is  perpendicular  to  the  orbital 
plane,  so  out-of-plane  thrusts  are  not  possible.  We  further  assume  the  chase  vehicle  is 
already  in  the  same  plane  as  the  target.  With  these  assumptions  we  will  develop  two 
different  solutions  to  the  Clohessy- Wiltshire  equations  and  evaluate  the  trajectories  and 
power  requirements  of  each  solution. 

5.1  Solutions  to  the  Clohessy~WUtshire  Equations 

In  section  2.5  we  introduced  the  Clohessy- Wiltshire  equations,  Ekiuations  (25)  and 

(26): 


X  —  2ny  —  3n*i  =  fr 
y  +  2nx  =  /» 

Our  approach  to  solving  the  forced  Clohessy- Wiltshire  equations  is  to  specify  a  de¬ 
sired  Xy  y,  i,  y,  X,  and  y,  then  use  Equations  (25)  and  (26)  to  determine  the  necessary  fr 
and  /«. 
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5.1.1  Decaying  Solution.  One  solution  is  to  require  x  and  y  to  behave  as  critically 
damped  oscillators,  i.e. 

x  =  e^'(h  +  rt)  (75) 

Then  the  differential  equation  of  the  system  is  necessarily  of  the  form 

i  +  ttii  +  o*x  =  0  (76) 

Write  the  characteristic  equation  and  solve  it  for  critical  damping: 

p*  +  oip  +  a,  =  0 


-Oi  ±  \/oJ  -  4a, 

^ - 2 - 

Setting  the  radical  equal  to  zero  establishes  critical  damping: 

oj  -  4a,  =  0 

1  2 
o,  =  -o? 

P-  -  2  2 

Next  we  return  to  Equation  (77)  and  its  derivative: 

X  =  c'*  (6  +  ci) 

X  =  (pb  +  pet  +  c)  (77) 

We  can  evaluate  Elquations  (77)  and  (79)  at  t  =  0  and  solve  for  fr  and  c  in  terms  of  the 
initial  conditions  x,  and  x,.  We  find  6  =  x,  and  c  s  x,  -  px,.  Thus  we  can  re-write 
Equations  (77)  and  (79)  as 


X  (t)  =  e"  (x,  -i-  (x,  -  px,)  t) 
x(0  =  c'‘(x,-»-p(x, -px,)0 
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Similarly  for  y: 


y(<)  =  e'‘(y.  +  (y,-py,)0  (80) 

y(0  =  «'‘(y.  +  p(y, -w.)0  (81) 


To  find  fr  and  /«  which  will  provide  the  desired  performance,  we  set  Equations  (25) 
and  (26)  equal  to  Equation  (78)  and  solve  for  f,  and  /»: 

X  -  2ny  -  3n’*  -  fr  =  i  +  Oii  +  c,x  =  0 

fr  =  -2ny  -  Oi*  -  (a,  +  i»*)  x  (82) 

y  +  2nx  -  /#  =  y  +  Oiy  +  a,y  =  0 

ft  =  2nx  -  Oiy  -  a,y  (83) 

where  x,  x,  y,  and  y  are  given  by  Equations  (80-83)  respectively.  The  value  0|  will  be 
chosen  depending  on  the  desired  system  performance. 
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5.1.2  Polynomial  Solution.  Another  approach  is  to  specify  Af  boundary  condi¬ 
tions  and  define  functions  x  and  y  as  (AT  —  l)th  order  polynomials  in  time.  For  instance, 
let 


X 

(84) 

X 

=  6, 2c,t 3d,t* 

(85) 

X 

=  2c,-»-6d,l 

(86) 

with  boundary  conditions 


x(0)  =  x,  x(f/)  =  0 
x(0)  =  x.  x(l/)  =  0 
*(</)  =  0 


If  we  evaluate  Equations  (86-88)  at  the  boundary  conditions  we  find 

a,  =  X, 

-3x,  6x, 

3x,  8x, 

T  T 

io  3x, 

In  a  like  manner 

y  =  a, b,t  + 
y  =  6,  -I-  2c,t  -1-  3d,P 
y  =  2c,  -f  6dft 
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(87) 

(88) 
(89) 


y(o)  =  y.  y(</)  =  o 
y(o)  =  y.  y(</)  =  o 
y(</)  =  o 


«t 

6, 

c. 


y. 


y. 

-3y. 

</ 


y»  3y. 
</  »/ 


Then  /r  and  /$  are  given  by  substituting  Ek|uations  (86-91)  into  Ex^uations  (25)  and 

(26). 


5.2  Methodology  of  the  Rendezvous  &  Docking  Study 

The  methodology  of  this  study  differs  slightly  from  chapters  3  and  4  due  to  the 
simplifications  we  introduced  at  the  beginning  of  the  chapter.  The  first  implication  is 
fl  =  {0  0  Bs}^  =  constant.  Then 


Solve  Equation  (92)  for  £  and  find 


(90) 
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This  simplifies  the  power  expression  also.  In  addition  to  the  lack  of  ks,  Br,  and  8$, 


we  can  simplify  the  relative  velocity  equation  to 


£r«j  =  (n-u^)rc0 


Thus  Equation  (47)  becomes 

P,  =  -«r  (VrtliBa)  +  4— («,*  +  K4^) 
m* 


As  in  chapter  3,  a  =  6, 578  km  for  the  reference  orbit,  which  has  a  period  T  =  5, 309*. 
We  choose  one  revolution  as  the  time  duration  in  which  to  perform  the  maneuver,  and  use 
tf  =  5,309*  in  evaluating  the  constants  of  the  polynomial  solution. 

An  unlimited  combination  of  initial  conditions  and  time  durations  can  be  studied. 
The  initial  conditions  we  studied  are  those  when  the  chase  vehicle  is  in  a  circular  orbit.  In 
that  case,  x,  =  0  for  any  x,  and  y,.  To  find  y,  we  revisit  Equation  (28): 

x(t)  =  -  f  ~  +  3x,^  cos  n<  +  “  sin  nt  +  4x,  +  — 

\  n  /  n  n 


If  X,  is  spedfied,  then  we  find  y,  by  setting  the  co^cient  of  the  cos  nt  term  equal  to  zero 
and  solving  for  in  terms  of  x,: 


n 


y« 


0 


SZgTl 

~~2~ 


oi  =  3n  is  used  in  the  decaying  solution  to  ensure  the  positions,  velocities,  and 
electrodynamic  accelerations  sufficiently  decay  within  tf.  oj  =  n  and  oi  =  2n  were  also 
considered,  but  neither  value  produced  x(t/)  and  y(t/)  sufficiently  close  to  zero  to  allow 
a  valid  comparison  between  the  decaying  solution  and  the  polynomial  solution.  For  ex¬ 
ample,  the  response  at  tj  of  the  decaying  solution  when  oj  =  2n  for  initial  conditions 
{xoiVo)  —  (1«1)  km  is  x{if)  =:  14  m  and  y(t/)  =  — 4  m.  In  contrast,  x{tf)  =  Tfim  and 
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Table  5.  Rendezvous  Power  Re>ulti>  {mlm^  —  10) 


Case 

Solution 

(km) 

y. 

(km) 

€ 

(J/kgxl0<) 

peak  P, 
(W/kg) 

P, 

(W/kg) 

1 

Decaying 

1 

0 

6.04 

234 

11.5 

2 

Polynomial 

8.65 

38.3 

16.5 

3 

Decaying 

1 

1 

1.63 

98.1 

3.10 

4 

Polynomial 

4.66 

22.1 

8.86 

5 

Decaying 

0 

1 

1.47 

14.5 

2.80 

6 

Polynomial 

.597 

2.68 

1.14 

7 

Decaying 

-1 

1 

14.2 

306 

27.1 

8 

Polynomial 

14.7 

61.4 

28.0. 

9 

Decaying 

-1 

0 

6.9o 

142 

13.1 

10 

Polynomial 

9.52 

40.1 

18.1 

11 

Decaying 

-1 

.1 

2.49 

52.1 

4.74 

12 

Polynomial 

5.52 

23.7 

10.5 

13 

Decaying 

0 

-1 

1.47 

60.6 

2.80 

14 

Polynomial 

.603 

3.81 

1.15 

15 

Decaying 

1 

-1 

13.4 

444 

25.5 

16 

Polynomial 

13.9 

59.5 

26.4 

y(tj)  =  -9/im  for  the  polynomial  solution  at  the  same  initial  conditions.  If  oj  =  3n,  then 
z  (ty )  s  0.9  m  and  y  (ty )  =  0.1  m,  which  we  deem  sufficiently  close  to  zero  to  allow  a  valid 
comparison  between  the  two  solutions. 


5.3  Results 

Both  solutions  were  evalu2,ted  at  eight  different  starting  positions  on  a  1  km  x  1  km 
square  centered  on  the  x  —  y  origin.  The  results  of  each  case  au-e  summarized  in  Tables  5 
and  6. 

The  decaying  solution  has  the  lowest  energy  and  P,  for  all  the  trajectories  except  the 
two  which  start  on  the  x  axis.  However,  in  every  case  the  decaying  solution  has  a  higher 
peak  P,  than  the  polynomial  solution.  If  peak  power  is  a  concern,  then  the  polynomial 
solution  appears  best,  but  the  tradeoff  is  sometimes  twice  times  as  much  total  energy  and 
P,,  depending  on  the  initial  conditions. 
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Table  6.  Rendezvous  Final  Position  and  Velocity  Results 


Case 

Solution 

*/ 

(m) 

Vl 

(m) 

(m/sec) 

»/ 

(m/sec) 

1 

Decaying 

.9254 

-.8357 

-.0015 

.0013 

2 

Polynomial 

7e-6 

-2e-5 

-3e-6 

7e-6 

3 

Decaying 

.9254 

.0897 

-.0015 

-.0002 

4 

Polynomial 

7e-6 

-9e-6 

-3e-6 

4e-6 

5 

Decaying 

0 

.9254 

0 

-.0015 

6 

Polynomial 

0 

7e-6 

0 

-3e-6 

7 

Decaying 

-.9254 

1.7612 

.0015 

-.0028 

8 

Polynomial 

-7e-6 

23e-6 

3e-6 

-lle-6 

9 

Decaying 

-.9254 

.8357 

.0015 

-.0013 

10 

Polynomial 

-7e-6 

16e-6 

3e-6 

-8e-6 

11 

Decaying 

-.9254 

-.0897 

.0015 

.0002 

12 

Polynomial 

-7e-6 

9e-6 

3e-6 

-4e-6 

13 

Decaying 

-9254 

.8357 

.0015 

-.0013 

14 

Polynomial 

0 

-7e-6 

0 

3e-6 

15 

Decaying 

.9254 

-1.7612 

-.0015 

.0028 

16 

Polynomial 

7e-6 

-23e-6 

-3e-6 

lle-6 

Each  solution  has  the  desirable  result  of  allowing  a  soft  dock.  The  velocities  and 
elec*  f-ynamic  propulsion  accelerations  all  go  to  zero  as  the  chase  vehicle  approaches  the 
target  in  every  case. 

Figures  37-40  show  the  trajectory,  velocity,  acceleration,  and  power  profiles  of  each 
solution  for  cases  1-4.  Note  the  peak  P,  occurs  at  the  beginning  of  the  thrust  period  for 
the  decaying  solution,  and  then  quickly  decays,  as  we  expect.  Both  solutions  cause  large 
power  discontinuities  at  In  each  case,  the  velocities  and  electrodynamic  propulsion 
accelerations  approach  zero  as  the  chase  vehicle  approaches  the  target.  Figures  37-40  are 
also  representative  of  cases  7-12,  15,  and  16. 

Cases  1-4  and  15-16  require  less  energy  than  their  counterparts,  cases  9-12  and  7- 
8.  This  is  to  be  expected,  since  cases  7-12  are  converting  electrical  energy  into  potential 
energy,  while  cases  1-4  and  15-16  are  trading  potential  energy  for  electrical  energy.  In 
cases  1-4  and  15-16  we  also  see  parts  of  the  trajectory  where  the  total  power  is  negative, 
indicating  the  ability  to  regenerate  power.  There  are  even  regions  in  the  first  quadrant 
{xo  >  0,  po  >  0)  where  the  total  energy  is  zero,  i.e.  the  entire  maneuver  can  be  done 
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Time  (sec) 


Figure  38.  Case  2:  P( 
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Specific  Power  (W/kg) 


Time  (sec) 


ial  Solution 


(oesAu)  A~a  X  , 


_ Table  7.  Adc 

Utional 

Rendezvous  Po 

.s 

Case 

Solution 

m/lTle 

t 

(J/kgxlO") 

peak  P, 
(W/kg) 

7 

8 

Decaying 

Polynomial 

6 

6 

8.72 

9.00 

156 

37.6 

■ 

7 

8 

Decaying 

Polynomial 

3 

3 

1 

8.7 

8.98 

with  no  net  energy  expenditure.  This  is  consistent  with  what  is  already  well  known  about 
electrodynamic  propulsion:  power  generation  creates  a  drag  which  lowers  altitude  (12:126). 
However,  if  a  mission  designer  accepts  the  constraints  on  initial  conditions  which  allow  the 
maneuvers  to  be  done  with  no  net  energy,  then  electrodynamic  propulsion  has  near  term 
application. 

The  worst  case  initial  conditions,  cases  7  and  8,  were  chosen  for  further  evaluation 
at  different  conductor  mass  rations.  These  results  are  summarized  in  Table  7.  These 
results  are  consistent  with  the  conclusions  in  Chapter  3:  the  /c?  terms  dominate  the  power 
equation,  so  if  larger  conductors  are  used  less  power  is  required.  This  worst  case  initial 
condition  is  still  beyond  the  4  W/kg  benchmark,  even  with  the  ‘aU  conductor’  ratio  of 
m/nif  —  3. 

Figures  41  and  42  show  how  either  solution  can  be  used  to  perform  a  V-BAR  ap¬ 
proach,  where  the  chase  vehicle  approaches  the  target  along  v.  fr  and  f$  are  controlled  to 
ensure  in-track  movement  without  any  radial  motion.  Both  trajectories  are  essentially  the 
same,  but  the  polynomial  solution  uses  less  energy. 

A  drawback  of  each  solution  is  the  power  discontinuity  at  t,.  This  can  be  solved  with 
a  new  fifth  order  polynomial  solution  and  the  additional  boundary  conditions 


z  (0)  =  z 


y(o)  =  y* 


Also,  even  though  Equations  (25)  and  (26)  are  coupled  in  z  and  y,  neither  the 
decaying  nor  polynomial  solutions  are  coupled.  If  y  can  be  specified  as  a  function  of  z  (or 
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Figure  42.  Case  13:  V-BAR  Approach  with  Decaying  Solution 
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vice  versa),  then  the  angle  of  approach  can  be  controlled.  Such  a  coupled  solution  might 
also  take  advantage  of  the  system’s  kinetic  and  potential  energy  in  such  a  way  as  to  lower 
the  electrical  energy  requirement.  Considerable  effort  was  expended  searching  for  such  a 
solution  without  success. 

Additional  study  should  be  performed  to  re-validate  these  results  without  the  sim¬ 
plifying  assumptions  used  here.  Also,  the  initial  conditions  were  very  specific,  and  are 
not  likely  to  be  found  in  actual  operations.  Various  initial  conditions,  especially  a  non¬ 
circular  initial  chase  orbit,  should  be  studied  to  ensure  the  conclusions  still  hold.  Out  of 
plane  motion  was  not  studied,  and  Lawrence  has  already  shown  undesired  out  of  plane 
motion  results  even  when  2,  =  0,  i,  =  0,  and  I,  =  0  (8:E-15,  E-17).  Future  work  should 
re-incorporate  this  fact. 
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VI.  Conclusions  &  Recommendations 

In  this  final  chapter  we  summarize  the  major  results  of  the  thesis  work  and  make 
recommendations  for  further  research. 

The  two  important  results  of  the  theoretical  development  were  the  definition  of  the 
vector  £  and  the  subsequent  derivation  of  the  specific  power  equation.  Together  they  allow 
the  analysis  to  proceed  without  having  to  specify  a  detailed  vehicle  configuration;  only 
designation  of  the  conductor  mass  to  total  vehicle  mass  ratio  and  selection  of  a  conductor 
metal  need  be  made.  The  specific  power  equation  clearly  shows  how  these  choices  affect 
the  power  required.  The  total  energy,  found  by  integrating  the  power  equation,  is  then 
useful  in  comparing  different  scenarios.  If  the  conductor  mass  to  total  vehicle  mass  ratio 
is  considered  analogous  to  a  classical  propulsion  system’s  propellant  mass  ratio,  then  the 
two  parameters  may  be  used  as  a  basis  for  comparison  between  electrodynamic  propulsion 
and  classical  propulsion  systems,  as  we  did  in  chapter  4.  In  the  future,  design  engineers 
can  use  specified  values  of  k  and  P,  to  perform  trade>off  studies  of  current,  conductor 
length,  or  number  of  conductor  turns  to  determine  the  optimum  configuration  for  a  specific 
application. 

Electrodynamic  propulsion  appears  to  be  a  feasible  near  term  alternative  to  classical 
propulsion  systems  for  small  orbital  plane  changes.  Plane  changes  are  best  done  at  rela¬ 
tively  high  inclinations  due  to  the  nature  of  the  magnetic  field.  Constraining  semi-major 
axis  change  may  expand  the  plane  change  rate  envelope. 

Electrodynamic  propulsion  may  also  be  used  for  stationkeeping  the  argument  of 
perigee  of  a  true  polar  Molniya  satellite.  The  Molniya  stationkeeping  strategy  developed 
here  is  very  sensitive  to  the  right  ascension  of  the  ascending  node.  A  candidate  for  future 
research  is  to  constrain  the  change  in  semi-major  axis  to  eliminate  this  sensitivity. 

For  both  plane  change  and  Molniya  stationkeeping,  development  of  forcing  function 
algorithms  incorporating  a  priori  knowledge  of  the  magnetic  field  and  application  of  op¬ 
timization  techniques  to  reduce  power  consumption  or  otherwise  improve  performance  is 
recommended. 
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Electrodynamic  propulsion  is  applicable  to  rendezvous  and  docking,  although  the 
relatively  high  power  requirements  limit  near-term  implementation  to  the  V-BAR  regime. 
If  spacecraft  specific  power  supplies  can  improve  to  approximately  25  W/kg  then  all  the 
rendezvous  cases  studied  in  Chapter  5  become  feasible.  The  ability  to  perform  a  soft 
dock  has  been  demonstrated,  which  is  advantageous  compared  to  a  two  burn  rendezvous 
and  docking  maneuver.  We  recommend  these  conclusions  about  rendezvous  aind  docking  be 
reconfirmed  after  incorporation  of  the  tilt  in  the  Earth’s  magnetic  field  into  the  geomagnetic 
field  model  and  consideration  of  out  of  plane  motion,  as  well  as  investigation  of  various 
other  initial  conditions.  Further  work  is  also  warranted  in  developing  solutions  to  the 
forced  Clohessy- Wiltshire  equations.  Since  the  Clohessy- Wiltshire  equations  are  coupled 
in  X  and  y,  we  recommend  investigating  a  solution  which  also  couples  x  and  y,  and  predict 
it  will  use  less  energy  than  our  uncoupled  solutions.  It  will  have  the  additional  benefit  of 
allowing  the  approach  angle  to  be  controlled. 

We  have  several  general  recommendations  for  future  research.  First,  the  feasibility 
of  power  supplies  which  can  actually  deliver  the  peak  specific  power  and  total  energy  in 
the  ranges  discussed  in  this  thesis  should  be  evaluated.  Second,  we  reiterate  Lawrence's 
recommendation  that  future  studies  be  conducted  with  a  more  accurate  geomagnetic  field 
model.  Third,  electromagnetic  compatibility  is  a  concern  for  the  integration  of  electro- 
dynamic  propulsion  into  any  vehicle,  and  should  be  studied.  A  particular  concern  is  the 
effect  of  the  localized  magnetic  field  caused  by  the  current  in  the  closed  circuit  conductor 
on  the  geomagnetic  field,  especially  if  the  vehicle  is  using  magnetometers  to  measure  the 
magnetic  field. 

The  most  feasible  near  term  realization  of  electrodynamic  propulsion  may  be  a  com¬ 
bination  of  chemical  and  electrodynamic  thrusters  applied  to  rendezvous  and  docking. 
The  chemical  thrusters  could  be  used  to  reduce  the  range  to  a  point  where  electrodynamic 
propulsion  can  take  over  and  perform  the  terminal  docking  maneuver  within  reasonable 
power  limits.  Perhaps  someday  electrodynamic  thrusters  can  even  be  incorporated  into 
manned  maneuvering  units,  allowing  astronauts  to  work  in  close  proximity  without  concern 
for  thruster  plume  impingement  on  each  other. 


Appendix  A.  File  Listings 

This  appendix  contains  the  Matlab  script  and  function  files  used  to  generate  the 
data  in  Chapters  3-5 

The  script  inc.m  was  used  for  the  inclination  change  study. 

%  variable  naaea  in  square  brackets  are  the  names  of  the  variables 

X  used  in  the  aain  body  of  tbe  thesis. 

tic 

V - - - —Get  started  — - — — — - 


i»i; 

n>1.18339e-3; 
alfao«0. : 


alf adot-7 . 2924d20S6e-S ; 

r>6S78000; 

Nag«8.05el5; 

fUAN-0; 

a»r; 
e«0. ; 
s»0. ; 

Mo»0; 

rhor»2.8e-8; 
rhoc»2700; 
mass*l ; 

Mu«39860i2el4; 

BiBC-10; 

idotbar«2 . 0257e-7 ; I*pi/2 ; 


Xidotbar>-2 . 0257e-07 ; I*pi/4 ; 
Xidotbar«2 .0257e-8 ; I«pi/2 ; 
Xidotbar*-2 . 0267e-08 ; I»pi/4 ; 


X  Counter  for  natrices  at  and 
X  Mean  action  (rad/sec) 

X  Inertial  longitude  of  Greensich 
Xaeridian  at  beginning  of  tine  span 
X  (rad)  Ctheta.g_0] 

X  Rotational  rate  of  Earth  (rad/sec) 
X  Cfflaega_o+] 

X  Initial  radius  (a) 

X  Magnetic  aoaent  [script  N] 

X  (kg/a*3-coul-sec) 

X  Right  ascension  of  ascending  node 
X  [cap  Oaega]  (rad) 

X  Saai-aajor  axis  (a) 

X  Eccentricity 
X  Arguaent  of  perigee  (rad) 

X  [little  oaega] 

X  Mean  anoaaly  at  epoch  (rad)  [H.o] 

X  Conductor  resistivity  (oha-a) 

X  [rho.R] 

X  Conductor  aass  density  (kg/a‘‘3> 

X  [rho.c] 

X  Total  vehicle  aass  (kg)  [a] 

X  Gravitational  paraaeter  of  the 
X  Earth  (a~2/s‘‘2)  [au] 

X  Conductor  aass  ratio  [a/a.c] 

X  idotbar  is  desired  average 
X  inclination  change  rate  (rad/sec) 

X  [bar  over  dot  over  I] 

X  I  is  inclination  (rad) 

X  1  deg/day 
X  I»45  deg,  -1  deg/day 
X  .1  deg/day,  I«90 
X  -.1  deg/day,  I»4S 
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pl«0;  t  Equinoctial  olMiont  [h] 

p2«0;  X  Equinoctial  alaaant  [k] 

M>0;  X  Naan  anoaaly  (rad) 

▼■0;  X  Arguaant  of  latituda  (rad) [u] 

- Initializa  na^atic  fiald  rotation  aatriz - 


phio>pi/2-78 . 3*pi/180 ; 
laao*69*pi/180 : 

RiaagBCcos(laBo)*cos(phio)  -ain(l«Bo)*coa(phio)  -sin(phio) 
sindano)  coa(lano)  0 

coa(lano)*8in(phio)  -8in(laao)*8in(phio)  coB(phio)]; 

X  Nota:  taag  rotataa  froa  Earth’s  cartasian  fraaa  to  aagnatic  cartasian 
Xfraaa  Raag^CR'gb] 


-  i,  tha  main  loop - 

for  t«to:t8:tf  X  t  is  tha  tiaa  (sac) 

X  ts  is  tha  tiaa  stap  (sac) 

X  tf  is  tha  final  tiaa  (sac) 

p«a*(l-a‘2);  X  Saai-latus  ractua  (a) 

Xif  raa(i,10)"«0  X  This  puts  out  a  aassaga  to  lat  you 

XCv* 180/pi  «* 180/pi]  X  know  vhara  tha  prograa  is  if  you 

Xand  X  «unt  it  to 

r»p/(l+a*co8(v)) ; 

alfaaalfao-t^alfadotot;  X  Rotatas  tha  Earth 


X - Croats  rot  aatriz  froa  orbital  fraaa  to  cart  inart 

rll«cos(RAAN)*cos(T)-sin(RAAI)*co8(I)asin(T) ; 
r 12*-cos (RAAN) *sin( v) -sin(RAAI) *cos (I) *cos(t) ; 
rl3»sin(RAAN)*8in(I) ; 

r21«8in(RAAN)*co8(T)'»co8(RAAN)*co8(l)*sin(T) ; 
r22*-sin(RAAN)*8in(T)^cos(RAAN)*cos(I)*cos(v) ; 
r23*~sin(I)*co8(RAAN) ; 

r31«8in(v)*8in(I) ; 
r32»8in(I)*co8(v) ; 
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r33«co«(I) : 


Ro2ci-Crll  rl2  rl3  X  CR'iaJ 


r21  r22  r23 
r31  r32  r33] ; 

X - Rotat*  [r  0  0]  in  sparical  coord  into  [z  j  z]  inartial 

zyzi«Ro2cl* [r  0  0] * : 

% - Rotata  xyzi  into  xyzE - 


Rci2cE«rot(0,alfa-alfao) :  X  Rci2cE  rotatas  froa  cartaaian 

X  inartial  fraaa  to  Earth’s 
X  rotating  cartaaian  fraaa  [R'gi] 

xyzE*Rci2cE*xyzi ; 


^ - ... —  rotata  xyzE  into  aagnatic  cartaaian  fraaa 

xyzaaRaagaxyzi ; 

- Coaputa  lat  and  long  in  aag  fraaa 


phia«atan2(xyza(3) ,  (xyza(l)'“2+xyza(2)*2)'.S) ; 
laBa*atan2(xyzB(2) ,xyza(l)) ; 


jj - - - - - Coaputa  aagnatic  fiald  coaponanta - - — 

z«2*Mag*ain(phim)/r*3;  X  Radial  coaponant  (kg/coul-aac) 

h«-Mag*coa(pbiB)/r‘3;  X  Tranavaraa  coaq>onant  (toaard 

X  north  pola)  (kg/coul>asc) 


X  -Coaputa  rot  aatrix  to  go  froa  apharical  aagnatic  to  cartaaian  aagnatic 
Raa2ca>rot(phia,laBa) :  X  CR'ba] 


X - Rotata  Cz  0  h]  into  cartaaiaa  aagnatic  fraaa 

Bca«Rsa2ca* [z  Oh]’;  X  (kg/coul-aac) 

X  -  Rotata  Bca  into  Earth’s  rotating  cartaaian  fraaa 

BcE^Raag’aBca: 
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X 


Rot At*  BcE  into  cartosian  inart ial  frana 


Bci«Rci2cE'«BcE;  X  (kg/cottl>sac) 

*yzI-Rci2cE**xyzE;  X  Dabug  coda:  zyzl«zyzi  if  all 

X  rotation*  ara  corract 

- Rotata  Bci  into  orbital  frana - 

B«Ro2ci'*Bci:  X  (kg/coul-sac) 

rtp«Ro2ci’*zy2l;  X  Dabug  coda  -  rtp»Cr  0  0]  if  all 

X  rotation*  ara  corract; 


- Coq>uta  kappa  vactor - 

b«B/nozB(B) ; 

d«Cb(2)  -b(l)  0]/*qrt(b(2)‘2+b(l)‘2);  X  k^ppa  unit  vactor  [kappa  hat] 


a3*idotbar*2*n*r*coa(T) ; 

iln»a3/*qrt(B(l)*2*B(2)“2) ; 
ILM«ilB*d; 

if  ■az(ab*(ILM}}>200 
ILM-200«ILM/Baz(ab*(ILN)) ; 
and 


X  Out  of  plana  accalaration 
X  (■/*ac“2)  [f.3] 

X  Scalar  il/n  (A-n/kg)  [kappa] 

X  IL/n  vactor  [kappa  vactor] 

X  Sorry,  no  infinite  kappa 


X  - -  Power  raquiranant*  - * - 

vnagE«[-alfadot*zyzE(2)  alfadot*zyzE(l)  0]';  X  B  field  velocity  [v.B] 

X  in  the  g  frana 


vnagi>Rci2cE  *  vvaagE ; 

vnago*Ro2ci ’ *vnagi ; 

vral« [r*a**in(v) *n/ ( l*a*co* (v) )  r*n 


X  v_B  in  the  i  frana 
X  v.B  in  the  a  frana 
0]’-vnago;  X  v.ral 


Pzi— (vral (2) *B (3) -vral (3)*B(2) ) ; 
Pyi»-(vral(3)*B(l)-vral(l)*B(3)) ; 
Pzi— (vral  ( 1 )  *B  (2) -vral  (2)  vB  ( I ) ) ; 
Pzr*4*rhor*rhoc*nnc* ILN ( 1 ) “ 2 ; 
Pyr«4*rhor*rhoc*nncvILM(2) “2 ; 
Pzr*4*rhor*rhoc*aac*ILM(3) “2 ; 
Pwri»Pxi+Pyi+Pzi ; 


X  iVi  radial 
X  iVi  in- track 
X  iVi  out  of  plana 
X  i“2R  radial 
X  i*2R  in-track 
X  i*2R  out  of  plana 
X  total  iVi 
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P»rr«Par+Pyr^Pir ; 


X  total  1*3R 


•uaP-CPwrl^Pwrr) :  X  total  [P.a] 

^  .... — - - - -  Croat  product  for  forcaa  — — - - - 

Fr«ILN(2)*B(3)>ILN(3)«B(2):  X  radial  forca  (k(-B/Bac‘2) 

Pu"XLN(3)«B(l)-IUI(i)*B(3);  X  Xn^tracB  forca  (kg-a/aac‘2) 

F3«XLN(1)*B(2}-XLM(2)*B(1):  X  out-of-plaaa  forca  (kf-B/sac‘2) 

\  — - - Collact  output  — — — - — - - 

X  1  2>4  5  e  7-9  10  11  12 

data(i.:)«Ctal80/pi  zysi*  phl*lB0/pl  laB*180/pl  B*  Fr  Fu  F3] : 

X  1  2-4  5-7 

datab(l,:)«ClaHi«180/pi  Boi*  Bel*]; 

X  12-45  6 

datia(l.:)«Cr*180/pi  zya*  laa*180/pl  phlB*180/pl] ; 

XdtaatCi. 180/pi  zysi*  zyzX*];  X  Dabuf  coda 

XdtaatCi. :)•[▼*  180/pi  10  0  rtp*]:  X  Dabuf  coda 

X  1  2345-78 

dat(i.:)»CT«180/pi  Fr  Fu  F3  XUf  ilmii 

X  1  2  3  4  5  6  7 

ala(i.:)*CT*180/pl  a  a  X*180/pi  IUAN*180/pi  ata2(pl.p2)*  180/pi  N*180/pi]: 
A3(i.l)»a3: 

P(i.:)»CPri  Pyi  Pai  Pxr  Pyr  Par  Pari  Pwrr  auaP]; 


X - Propogata  alaaata - 

lagraoca;  X  Subrout iaa  Lagranga  iatagratas 

X  Lagraaga’a  plaatary  aquatioaa 

X  Updata  atriz  iadaz 
aad 
toe 


Bubplot(2.2,l} 

plot(ala< : .  1)  .P( ; ,  1) .  ’ :  ’  ••Imi ;  .1)  .P( ;  ,2) ,  •  .ala( ; .  1)  ,P( :  .3) ,  •  —  * , 

X  ala(:  .1),P(:  ,4)  ,*-*)  X  Thia  liaa  ia  a  coatiauatioa  of  tha 

X  prarioua  oaa 
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jlab«l( 'Specific  po««r  (H/kg)') 

subplot (2, 2, 2) 
plotCslsaC : . 1) .slsaC : .4)) 
ylsbsK  *  Inclination*) 

subplot (2, 2, 3) 

plot(sls«( ; .1) ,dat( : .2) , *y* .sla*(: .1) ,dat(: .3) , *r* .s1ob( : .i) .dat( : ,4) . *g* , 

X  s1ob(: ,1) ,A3, *w*)  X  Ibis  lino  is  a  continuation  of  tho 

X  provious  ono. 

Xplot(dat ( : .1) ,dat( : .4) .dat( : .1) .13) 
ylaboK 'Accol oration  (■/ 000*2) * ) 

subplot (2. 2. 4) 

plot(olosi(:,i).dat(:.S).olsa(:.i).dat(:.6).*r'.olsa(:.l).dat(:.7),'g*) 

ylabolCIL/n*) 


■■sizoCP.l);  X  lu^or  of  rows  in  P 

onorgy«(.6*(P(l,4)+P(*.4))+su»(P(2:«-l,4)))*ts  X  Trapozoidal  intogration 
nazP«Baz(P) 

avgP»onorgy/(tf-to)  X  Cbar  osar  P.s]  (W/kg) 

da«a>6S78000  X  Chango  in  a  (■) 

I*I*180/pi 
disp(*bit  a  koy*) 
pauso 

subplot (2. 2,1) 
plot(oloa( : ,1) ,oloa( : ,5)) 
ylabol('RAAM’) 

subplot(2,2,2) 
plot(oloB(: ,1) .dobugC; ,4)) 
ylaboK’h*) 

subplot (2. 2, 3) 
plotColoaC : , 1) ,oloa( : ,2)) 
ylabol(*saBi-Bajor  axis’) 

subplot (2, 2, 4) 

plot(olaB(:.l),olsa(:,3)) 

ylabol('k') 
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The  script  reaa.a  was  used  for  fl  change  study. 

X  variable  names  in  square  brackets  are  the  names  of  the  variables 

X  used  in  the  main  body  of  the  thesis. 

tic 

jj  „ — - — _ — started  —————— — — — - 


i*l; 

n>l .  18339e-*3 : 
alfao>0. ; 


alf adot-7 . 2924620S6e-S ; 

r-6578000; 

Mag«8.0SeiS: 

RAAH-0; 

a»r; 

e«0. ; 

w»0. ; 

No«0; 

rhor«2.8e*8: 
rhoc»2700 ; 

mass*l ; 

Mu-3986012el4; 

nmcalO; 

RAARdotbar«2 . 02S7e-7 ; I«pi/2 ; 


XRAAIdotbar— 2 . 0257e-07 ; I-pi/4 ; 
XRAAVdotbar«2 . 0257e-8 ; I-pi/2 ; 
XRAARdotbar— 2 . 0257e-08 ; I-pi/4 ; 
pl-0; 
p2-0; 

M-0; 

v-0; 

X 


X  Counter  for  matrices  at  end 
X  Mean  motion  (rad/sec) 

X  Inertial  longitude  of  Greenwich 
Xmeridian  at  beginning  of  time  span 
X  (rad)  Ctheta.g.0] 

X  Rotational  rate  of  Earth  (rad/sec) 
X  [omega.o-»] 

X  Initial  radius  (m) 

X  Magnetic  moment  [script  M] 

X  (kg/m“3-coul-sec) 

X  Right  ascension  of  ascending  node 
X  [c^  Omega]  (rad) 

X  Sead-major  axis  (m) 

X  Eccentricity 
X  Argument  of  perigee  (rad) 

X  [little  omega] 

X  Mean  anomaly  at  epoch  (rad)  [N.o] 

X  Conductor  resistivity  (ohm-m) 

X  Crho_R] 

X  Conductor  mass  density  (kg/m‘3) 

X  [rho.c] 

X  Total  vehicle  mass  (kg)  [m] 

X  Gravitational  parameter  of  the 
X  Earth  (a‘2/s‘2)  [mu] 

X  Conductor  mass  ratio  [m/m_c] 

X  RAAldotbar  is  desired  average 
X  RAAM  change  rate  (rad/sec) 

X  [bar  over  dot  over  Omega] 

X  I  is  inclination  (rad) 

X  True  polar  sun  synchronous  case 
X  1-45  deg  sun  synchronous 
X  .1  deg/day,  1-90 
X  -.1  deg/day,  1-45 
X  Equinoctial  element  [h] 

X  Equinoctial  element  [k] 

X  Mean  anomaly  (rad) 

X  Argument  of  latitude  (rad)[u] 


Initialize  magnetic  field  rotation  matrix 
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phioT»i/2-78 . 3*pi/180 ; 
laBO«69*pi/180; 

RMg"Ccos(laBo)*cos(phio)  -■in(laao)*coa(phio)  -ainCphio) 
aindaan)  coa(laao)  0 

coa(laBo)*ain(phio)  -aln(laBo)*ain(phio)  coa(phio)] ; 

X  Iota:  Haag  rotataa  froa  Earth* a  cartaaian  fraaa  to  aagnatic  cartaaian 
Xfraaa  RaagaCR^gb] 


H -  wAin  loop - 

for  tato:ta:tf  X  t  ia  tha  tiaa  (aac) 

X  ta  ia  tha  tiaa  atap  (aac) 

X  tf  ia  tha  final  tiaa  (aac) 

paaad-a'Q) :  X  Saaidatua  ractua  (a) 

Xif  raB(i,10)»0  X  Thia  puta  out  a  aaaaaga  to  lot  you 

XCt*i80/pi  w* 180/pi]  X  )moa  nhara  tha  prograa  ia  if  you 

Xand  X  want  it  to 

r«p/(l*aacoo(T)); 

alf a"alf ao-^alf adot*t :  X  Rotataa  tha  Earth 


X - - —  Craata  rot  aatriz  froa  orbital  fraaa  to  cart  inart 

r 1 l*coa (RAAI) «co8 (▼) -ain(RAAI) *coa (I)*ain(T) ; 
rl2>-coa(RAAI)*8in(T)-ain(RAAI)*coa(I)*coa(r): 
rl3«-8in(RAAII)*8in(I) ; 

r21*ain(RAAH)*co8(T)-»co8(RAAR)*co8(I)«8in(T): 
r22*-8in(RAAM)  *8in(T)«co8  (RAAI)  acoad)  *co8  (▼) ; 
r23«-8in(I)*co8(RAAI) ; 

r31*8in(T)*8in(I) ; 
r32>8in(I)*co8(v) ; 
r33«co8(I); 

Ro2ci>[rll  ri2  rl3  X  CR'ia] 

r21  r22  r23 
r31  r32  r33] ; 
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X - Rotat*  [r  0  0]  in  iparicnl  coord  into  [x  j  s]  inortial 

z7zl«Ro2ci*Cr  00]’; 

j - Rotata  xyzi  into  xyzE - 


Rci2cE>rot(0,alfa>alfao) ;  X  Rci2c£  rotataa  froB  cartasian 

X  inartial  frasa  to  Earth’s 
X  rotating  cartasian  fraaa  CR‘gi] 


xyzE«Rci2cE*xyzi : 


X  — rotata  xyzE  into  nagnatic  cartasian  fraaa  — 
xyzB«Raag*xyzi ; 

- - - - Coaputa  lat  and  long  in  nag  fraaa  — - 

phiB*atan2(xyzB(3) , (xyzB(l)‘2'»xyzB(2)‘2)‘.5) ; 
lanBaatan2(xyzai(2)  .xyzaCD) ; 

j - - — - - Coi^uta  Mgnatic  fiald  co^>onants 


z«2*llag*sin(phiB)/r‘3:  X  Radial  coaponant  (kg/coul-sac) 

h«-Mag*cos(phiB)/r‘3;  X  Transvarsa  coaponant  (toward 

X  north  pola)  (kg/coul-sac) 

X  -C<»q>uta  rot  natrix  to  go  froB  spharical  aagnatic  to  cartasian  aagnatic- 
RsB2cB"rot(phiB,la«) :  X  [R~bB] 

X  — - - Rotata  Cz  0  h]  into  cartasiaB  aagnatic  frana - 

Bcb«Rsb2cb*  Cz  Oh]*;  X  (kg/coul-sac) 

X - Rotata  Bcb  into  Earth’s  rotating  cartasian  frasa - 


Bc£>Raag’*BcB; 

jj - Rotata  BcE  into  cartasian  inartial  fraaa 

BciaRci2cE’*BcE;  X  (kg/coul-sac) 


xyzI«Rci2cE’*xyzE:  X  Dabug  coda:  xyzl»xyzi  if  all 

X  rotations  ara  corract 
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X 


Rotat*  Bel  into  orbital  fri 


B>Ro2ci**Bci;  X  (kg/coul>sac) 

rtp>Ro2ci**zjzI;  X  Oabug  coda  -  rtp»[r  0  0]  if  all 

X  rotations  ara  corract; 

X - - - - - Coi^uta  kappa  vactor - 

b«B/non(B) : 


d-Cb(2)  -b(l)  0]/aqrt<b<2)“2'*b(l)*2); 
a3aRAABdotbar«2*n*a*sin(I)*sin(T) ; 


iln-a3/aqrt(B(l)“2+B(2)“2) ; 
ILN-iln*d: 

if  Baz(aba(ILM))>200 
ILM-200*IUf/aAz(aba(ILM)) : 


vaagi«Rci2cE '  *magE ; 

Tnago«Ro2ci  ’  *vnagi ; 

▼ral«Cr*a*8in(r)*n/(l«a*cos(T))  r*n  03 

Pxi— (▼ral(2)*B(3)-rral(3)*B(2)) ; 
Pyi«-(vral(3)*B(l)-rral(l)*B(3)); 
Pzi— (Tral(l)*B(2)-Tral<2)*B(l)); 
Pzr«4*rhor*rhoc*nmc*ILK ( 1} *2 ; 
P7r>4*rhor*rhoc*nac*ILN(2) “2 ; 
Pzz«4*rhor*rhoc*nBc*ILM(3} “2 ; 
Pwri»Pzi+Pyi+Pzi ; 

Pwrr»Pzr ♦Pyr+Pzr ; 

sunP«<Psri+P¥rr) ; 

jj - Cross 


X  kappa  unit  tactor  [kappa  bat] 

X  Out  of  plana  accalaration 
X  (■/sac‘2)  Cf.3] 

X  Scalar  il/m  (A-n/kg)  [kappa] 

X  ll/m  Tactor  [k^pa  ractor] 

X  Sorry,  no  infinita  kappa 


X  B  fiald  ralocity  [r.B] 
X  in  tha  f  fraaa 

X  r.B  in  tba  i  fraaa 

X  r.B  in  tha  a  fraaa 

-vaago;  X  r.ral 

X  iVi  radial 
X  iVi  in- track 
X  iVi  out  of  plana 
X  i‘2R  radial 
X  i*‘2R  in-track 
X  i"2R  out  of  plana 
X  total  iVi 
X  total  i*2R 

X  total  [P.s] 

product  for  forcas  — — - - - 


and 

X - — —  Povar  raquiraaants 

vaagE>[-alfadot*xyzE(2)  alfadot*zyzE(l)  0]»; 
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rr«ILII(2)*B(3)-ILI!(3)*B(2): 
Fu-ILM(3)*B(1)-ILI((1)*B(3) ; 
F3>ILN(1)*B(2)>ILN(2)*B(1) : 


X  radial  forca  (kg-m/sac‘‘2) 

X  In-track  forca  (kg-B/aac‘2) 

X  out-of-plana  forca  (kg-B/aac*2) 

j - - - — — —  Collact  output  — — - 

X  1  2-4  5  6  7-9  10  11  12 

data(i,;)»[v*180/pi  xyzi*  phi*180/pi  laB*180/pl  B*  Fr  Fu  F3] ; 

X  1  2-4  S-7 

datab(i.:)*ClaBa*180/pi  Bea'  Bci']; 

X  12-45  6 

dataaCi, :}«Cv*180/pi  zjza*  laaa«180/pi  phla* 180/pi] ; 

XdtastCi, :)*Cv* 180/pi  zjzi*  zyzl*];  X  Dabug  coda 
XdtaatCi, :)•[▼*  180/pi  10  0  rtp*]:  X  Dabug  coda 

X  1  2345-78 

dat(i.:)«CT*180/pi  Fr  Fu  F3  ILM  ila] ; 

X  1  2  3  4  5  6  7 

alaa(i, :)«Cv*180/pi  a  a  I*180/pi  RiAM*180/pi  atan2(pl,p2)*180/pi  M*180/pi] ; 
A3(i,l)*a3; 

P(i,:)»CPxi  Pyi  Pzi  Pzr  Pyr  Pzr  Pwri  Parr  suaP]; 


- - - -  Propogata  alaaanta - 

lagranga;  X  Subrout ina  Lagranga  intagrataa 

X  Lagranga ’a  planatary  aquations 

i*!*!;  X  Updata  aatriz  indaz 

and 

toe 


subplot (2, 2,1) 

plot(alaa(:,l),P(:,l).';*,alaa(:,l),P(;,2),'-.',alaa(;,l).P(;,3),»— 

X  alaaC: .1),P(: ,4) X  This  lina  is  a  continuation  of  tha 

X  prarious  ona 

ylabaK’Spacific  powar  (H/kg)*} 

subplot (2, 2, 2) 
plot(alaa(: ,1) .alaaC: ,4)) 
ylabaK  ’  Inclination’ ) 
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subplot (2. 2, 3) 

plot(slM(:  .l),dst(:  .2).*7'  ,sl«i(: .1) .dstC:  .3)  . 'r'  ,s1sb( :  ,1) ,dat(:  ,4)  . , 

X  slsa(: ,1) ,A3, '«*)  X  This  lino  is  s  continuation  of  ths 

X  prsvious  ons 

Xplot(dat(: .1) .dat(: .4),dat(: .1).A3) 
ylabsK ' Accslsration  (B/tac*2) ' ) 

subplot (2. 2. 4) 

plot(sloB( : . 1) .date : .5) .slsa( : .1) .dat( : .«) . *r ' .slsaC : . i) .dat( : .7) . ) 
ylabsK'IL/n*) 


■•sizs(P,l):  X  lu^sr  of  rows  in  P 

snarg7a(.5*(P(l,4)-»P(B,4))-»suB(P(2:B*1.4)))*ts  X  Tr^azoidal  intsgration 
naxP«Baz(P) 

avgP>ansrg7/(tf-to)  X  [bar  osar  P.s]  (H/kg) 

da^a^dSTSDOO  X  Changa  in  a  (■) 

I>I* 180/pi 
dispC'hit  a  kay') 

pausa 

subplot (2, 2,1) 

plot(alaB(:,l),alaB(:.S)) 

ylabal(*IUAI'} 

subplot (2, 2, 2) 
plotCalaaC: ,1) .dabugC: ,4)) 
ylabaK'h*) 

subplot (2, 2, 3) 
plotCalaaC : , 1) .alaaC : ,2)) 
ylabolC'SMii -major  axis’) 

subplot (2, 2, 4) 
plot(alam(; ,l),alam(: ,3)) 
ylabal(’k’) 


91 


The  script  lagrangelr.a  is  called  by  inc.B  and  raan.a.  It  evaluates  Lagrange’s 
planetary  equations  and  integrates  the  orbital  dements  used  in  the  plane  change  studies. 

X  Thin  filn  i^;>lmnnts  Lagrangn’s  planetary  equatlona 
X  Square  brackets  denote  the  noaeclature  used  in  the  thesis  writeup 

adot«2*e*sin( v) / (n*sqrt ( l-e'’2) ) *Fr/Basa«2«aesqrt ( l-e*2}/ (n*r) eFu/aass ; 
Idot«  (  (r*cos  (s+w  )  )  /  (n*a*a*sqrt  ( l-e*2)  )  )  eFS/uss ; 

IUAMdot«(  (r*sin(w-»w)  )/ (n«a*a*sqrt(  i-e"2)*sin(I)  ) )  *F3/Bass ; 

dM«n-2*Fr/(a*n*Bass)-sin(theta)*F3/(tan(I)*a*n*Bass) ;  X  This  is  [ds/dt] 

X  pi  is  Ch] 

X  Next  line  is  first  two  texas  of  the  h  dot  eqtn 

pldotru>'‘Cos(theta)*Fr/(n*a*aass)'»(pl*r-*^(r'»a)*sin(theta))*Fu/(n*a‘2*aass) ; 

X  Next  line  is  the  last  tera  of  the  h  dot  eqtn 
pldot3a'p2**r*sin(theta)*F3/(n*a“2*tan(I)): 

pldotapldotru'*’pldot3:  X  This  is  all  three  terns 

X  p2  is  [k] 

X  Next  line  is  the  first  two  terns  of  the  k  dot  eqtn 

p2dotru«sin(theta)*Fr/  (n*a*Bass)<*'(  (r«a)*cos  (theta)  ♦p2*r)*Fu/  (n*a‘2*Bass) ; 

X  Next  line  is  the  third  tern  of  the  k  dot  eqtn 
p2dot 3>p 1 *r *s in (thet a) «F3/ (nea" 2*t an (I ) ) : 

p2dot«p2dotru-*^p2dot3 :  X  This  is  all  three  teras 

Xdn*- 1 . 5*sqrt (Mu/a‘S) *adot ; 

X  Next  lines  are  the  Euler  integration  of  eleaents 

pl*pl-»pldot*ts; 

p2>p2'*^p2dot*ts ; 

a«a-»adot*ts ; 

e*sqrt(pl“2+p2‘2) ; 

vadHets^w;  X  This  is  [s] 

I*I<»Idot*ts ; 

RAAN*iUAN4RAANdot*ts ; 

debugCi,  :)«Ci  t  reaCw*  180/pi, 360)  pi  p2  sqrt(pl‘'24p2*2)] ;  X  Debug  data 

X  collection 
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The  script  aolaiya.B  was  used  for  the  Moloiya  statioakeeping  study: 


tic 

X  Variable  nasas  in  equare  brackets  are  those  used  in  the  sain  body  of 
X  the  thesis. 

X  This  is  the  file  used  to  investigate  the  following  designer  orbit: 

X  inclination:  90  deg 
X  period:  86160  s  (*24  hr) 

X  rotation  of  arguaent  of  perigee:  0  d^day 

^ - Q^x  Started - 


i-1;  X 

n>2*pi/86160  ;  X 

alfao"0.;  X 

X 

X 

alfadot-7.2924620S6e-S:  X 

X 

Mag«8.0Sel5;  X 

X 

I-pi/2;  X 

aass«l;  X 

Mu«3.986012ei4:  X 

e-. 82501;  X 

v*1.5*pi:  X 

H»pi;  X 

ke— 2.67542e25;  X 

X 

X 

vdot«0;  X 

rhor«2.8e-8;  X 

rhocB2700:  X 

wdotbar»-3.8e-7  X 

X 

X 

X 

X 

RAANo-RAAH  X 


Counter  for  aatrices  at  end 
Hean  Motion  (rad/sec) 

Inertial  longitude  of  Greenwich 
■eridian  at  beginning  of  tiaia  span 
(rad)  [little  QBega_oplus_o] 
Rotational  rate  of  Earth  (rad/sec) 
[little  oaiega.oplus] 

Magnetic  aosent  (kg*a'‘3/coul-sec) 
[script  M] 

Inclination  (rad) 
mass  (kg)  M 

gravitational  paraaeter  (a‘3/s*2); 
Eccentricity 

Arguaent  of  perigee  (oaega)  [rad] 
Mean  anoaaly  (rad) 
term  for  Earth’s  oblateness  equal 
to  G  times  mass  of  the  Earth 
see  section  4.1  (units?) 
oaega  change  rate  (rad/sec) 
Conductor  resistivity  (oha-a) 
Conductor  aass  density  (kg/a‘3) 
Desired  «»ega  change  rate 
(rad/sec).  This  value  has  to  be 
adjusted  to  obtain  the  desired 
change  in  oaega.  The  noainal  value 
is  wdotbara-S.dOCe-T  rad/sec 
RAAIo  is  intitial  RAAM  (rad) 


^ - - - lait  field  rotation  aatrix 

phio-pi/2-78 . 3*pi/ 180 ; 
lano«69*pi/180; 
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RBag«[coa(laao)*cos(phio)  -■in(laBo)*coa(phlo)  -ainCphio) 
sinClaBo)  cosCluo)  0 

cos(laBo)*ain(phio)  -ain(laBo)«ain(phio}  coa(phio)]; 

X  lota:  Rug  rotataa  froa  Earth* a  cartaaian  fraaa  to  aagnatic  cartaaian 
X  fraaia.  Raag  ia  CR‘gb] 


j - - — _ — - —  ig  tha  main  loop  - - - - - - — - 

for  t«to:ta:tf  X  t  ia  tha  tiaa  (aac) 

X  ta  ia  tha  tiaa  atap  (aac) 

X  tf  ia  tha  final  tiaa  (aac) 

p«a«(l-a‘2):  X  Saai-latua  ractta  (■) 

Xif  raB(i,10)»0  X  Thia  puta  out  a  naaaaga  to  lat  you 

X[t« 180/pi  a* 180 /pi]  X  ahara  tha  prograai  ia  if  you  aant 

Xand  X  it  to 

EoaM;  X  Bagin  algoritha  to  aolaa  for  E 


E«kaplar(M.Eo,a) ; 
for  j«l:10 
if  aba(E-Eo)>1.0E-9 
Eo*E: 

E*kaplar(N,E,a) ; 

alaa  and 

and  X  End  algorithn  to  aolva  for  E 

▼■acoa((o-coa(E))/(a*coa(E)-l)) ;  X  Trua  anomaly  from  E 

if  ain(N)<0  X  Quadrant  chack 

and 


r«p/(l-»a*coa(T)) :  X  r  calculation 

▼aat*Ca*ain(v)*aqrt(Mu/p)  (l'»a*coa(T))*aqrt(Mu/p)  0];  X  aalocity  in  a  fraaa 

alfaaalfao-^alfadotat;  X  Rotataa  tha  Earth 

X - Craata  rot  Duitriz  from  orbital  frama  to  cart  inart - 

r 1 l«coa (RAAH) acoa (vav) -ain(RAAI) *coa (I)*ain(va«) ; 
rl2«'Coa(RAAR)*ain(vaT)-ain(RAAI)*coa(I)*coa(T-«>«) ; 
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rl3a-sin(RAAI)*iin(I) ; 


r21«sin(IUAN)*cos(v««)-»coB(RAAI)*cos(I)*ain(v««); 
r22«-8in(IlAAI)*ain(T'»«)<»cos(IUAI)«coa(I)*cos(v‘»>) ; 
r23«-sin(I)*cos(RAAN) ; 

r31»«in(T*w)*«in(I) ; 
r32«siii(I)*cos(v-»«) ; 
r33acoa(I) : 

Ro2ci*Crll  rl2  rl3 
r21  r22  r23 

r31  r32  r33] ;  X  CR‘ia] 

X - — ' —  Rotata  [r  0  0]  in  aparical  coord  into  [z  j  z]  inartial 

zyzi«Ro2ci* [r  0  0] ' ; 


- — - —  Rotata  xyzi  into  xyzE  ———————— - 

Rci2cE>rot(0,alfa'-alfao) :  X  Rci2cE  rotataa  froa  cartaaian 

X  inartial  fraaa  to  Earth’a  rotating 
X  cartaaian  Iraaa  £R‘gi] 


zyzEaRci2cE*zyzi ; 

% - —  coi^uta  lat  and  long  in  g  fraaia - 

phiE*atan2(zyzE(3) ,(zyzE(l)*2'»zyzE(2)'2)*.5) ;  X  latitude  (rad) 

laiaEBatan2(zyzE(2) .zyzE(l)) ;  X  longitude  (rad) 


1  - — -  rotata  zyzE  into  aagnatic  cartaaian  fraaa  - 

zyza«Raag*zyzi ; 

- - - -  Co^uta  lat  and  long  in  sag  frana 


phiasatan2(zyzm(3) ,  (zyzB(l)*2'»zya(2)'2)*.5) ; 
lajim«atan2(zyzm(2)  .zyzad)) ; 

% - Co^>uta  aagnatic  field  coaiponanta - 

z*2*Nag*ain(phim)/(r'‘3) ;  X  Radial  coaqponant  (kg/coul-aac) 
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h«”M«g*co«(phiM)/(r“3) ; 


t  Transvars*  coconut  (tovard 
X  north  pola)  (kg/coul-aac) 

X  Saa  Tasciona  for  dataila 


X-  Conputa  rot  aatriz  to  go  from  spharical  aagnatic  to  cartaaian  aagnatic 
R8ai2ca*rot(phim,lanm) ;  X  [R^bn] 


X - Rotata  [z  0  h]  into  cartaaian  nagnatic  frana 

Bcm>R8a2cB* [z  0  h] * ;  X  (kg/coul>sac) 

X  -  Rotata  Ben  into  Earth’s  rotating  cartaaian  frana 


BcE«Rnag’aBcn; 


X 


Rotata  BcE  into  cartaaian  inartial  frana 


Bci«Rci2cE'*BcE:  X  (kg/coul*sac) 

zyzI«Rci2cE**zyzE:  X  Dabug  coda:  zyzl>zyzi  if  all 

X  rotations  ara  corract 

^ - Rotata  Bci  into  orbital  frana - 

B»Ro2ci**Bci;  X  (kg/coul-aac) 

rtp«Ro2ci’*zyzI:  X  Dabug  coda  -*  rtp«[r  0  0]  if  all 

X  rotations  ara  corract; 

- Conputa  kappa  factor - - 

fr»(kaa(l-3*(sin(ii+v)*sin(I))“2))/(r“4) ;  X  radial  accalaration  (n/sac~2) 
fu»(ka*ain(2*(w+v))*(ain(I)“2))/(r“4);  X  in-track  accalaratoin  (n/sac“2) 
fl»ka*sin(w»T)*sin(2*I)/(r*4) ;  X  out  of  plana  accalaration  (n/s“2) 

sub.adoto  X  Calculatas  tha  affact  of  tha 

X  Earth’s  oblatanass 

^ - Co^>uta  IL/n  Tsetor - 

ILM-CO  0  0]’; 

if  cos(v)>0.1;  X  Linit  thrusting  to  this  rsgion 

al*8qrt(l-a‘2)/(n*a*a) ;  X  [k_l] 

bl«l+l/(l+a*co8(v)) ;  X  Ck.2] 
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X  iMt  t«o  lines  are  [kappa. that s] 

ILM(2)«-wdotbar/(-al*bl~2*sln(v)*2aB(3)/cos(v)-al*cos(T)*B(3)- 
al*(cos(T)*B(2)<»bl*sin(T)*B(l) )  “2/(coa(T)*B(3)  )  ) ; 

ILll(l)«bl*tan(v)*ILM(2) :  X  Ck^pa.r] 

ILM(3)*-ILM(2)*(cos(T)«B(2)-»bl*sln(T)*B(l))/(cos(v)*B(3)):  X  Ck^pa.3] 

ILNckal.01*ILM:  X  (Slack  to  anaura  function  has  boon 

if  noza(ILMck)-noni(IlJ()<0  X  niniaizad 

disp(*Waxning  -  non-niniatm  posar  at*) 

[i  t  TaiBO/pi] 
and 

if  ■az(abs(IIJI))>500  X  Liait  aagnituda  of  kappa  if 

ILN-SOO*ILM/aax(abs(ILN)):  X  nacassary 

and 

and 

X  — - - Croa-  product  for  forces  — — - 

F1«IUI(2)*B(3)-IL)((3)*B(2):  X  radial  force  (kg>a/sac‘2) 

F2>ILN(3)*B(1)-ILN(1)*B(3);  X  In>track  force  (kg*a/sac‘2) 

F3«ILN(1)*B(2)-ILN(2)«B(1):  X  out-of-plana  force  (kg-a/sac*2) 


Fr*Fl-*fr;  X  Cosines  affects  of  Earth 
Fu"F2afu:  X  oblatanass  and  thrusting 
Fl*F3afl;  X  fr,  fu,  and  fl  C(»a  froa 
XFu«2afu;  X  sub.vdoto 

X  — - - Power  raquiraaants  — — - — — - — — - 


TaagE*[-alfadot*z7zE(2}  alfadot*zyzE(l)  0]*;  X  B  field  velocity  [v.B] 

X  in  the  g  fraaa 

vaagi>Rci2cE’*vaagE;  X  ▼.B  in  the  i  fraaa 

vBagoBllo2ci*aTaagi:  X  w.B  in  the  a  fraaa 

vral«vsat ’ -vaago ;  X  w.ral 

X  Next  two  lines  are  the  power  due  to  the  induce  voltage 

Pi— <(vral(2)*B(3)-vral(3)*B(2))aILM(l)*(vral(3)*B(i)-vral(l)*B(3))*ILII(2) 

•»(vral(l)*B(2)>vral(2)*B(l))«ILM(3)); 

X  Iczt  3  lines  are  the  power  due  to  resistance  in  the  radial,  in-track, 

X  and  out  of  plane  conductors  respectively 
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Px»-(»r«l(2)*B(3)-vr«l(3)*B(2))*IIJI(l)+4*rhor*rhoc**ic*ILM(l)‘2; 
Py«-(Tr*l(3)*B(l)-»r*l(l)*B(3))*IlJI(2)+4*rhor*rhoc*«ic*ILII(2)‘2; 
Pz»-(»r«l(l)*B(2)-vr*l(2)*B(l))*ILII(3)+4*rhor*rhoc**ic*ILM(3)‘2; 
■u»P«(P*>Py+Pz) ; 


- -  Collect  output - 

X  1  2-4  5  6  7-9  10  11  12  13 

data(i,l:13)«CT*180/pi  zyzi'  phiE*180/pi  laBE*180/pi  B'  thata*180/pi  Fr  Fu  F3] ; 

X  1  2-4  6-7 

datab(i,l:7)a[laaun*180/pi  Bca’  Bci*]; 

X  12-45  6 

datam(i,l:6)aCv*180/pl  zyza*  laBi*180/pi  phia* 180/pi] ; 
dtast(i.l:7)aCT*180/pi  zyzi*  zyzI*];  X  Dabug  coda 
dta8t(i,l:7)«Cv*180/pi  100  rtp'/r];  X  Dabug  coda 

P(i.:)*CPi  Pz  Py  Pz  auaP] ; 

X  1  23  4  5  678 

datCi. 1: 8) -[▼* 180/pi  t  I*180/pi  iUAI*180/pi  (1.5*pi-«)* 180/pi  Fr  Fu  F3] ; 

X  12  3  456789  10  11 

aJigla8(i,l:4)aC  t  raiOO/pi  N*180/pi  E*180/pi]i 
kap(i.l:8)aC  t  talOO/pi  a  a  I  RAAI  «  Mo]; 

f orcaad.l; ID- [t  180/pi  fr  fu  fl  Fl  F2  F3  Fr  Fu  FI]  ; 

curraiit(i,l:5)-[t  v*180/pi  ILN’] ; 

Hdot(i,l:3)-C«dot  vdoto  «dot-»«doto] ; 

alanCi, :)>Ct  a  a  I*180/pi  RAAN* 180/pi] ;X  «*180/pi]:X  No* 180/pi] ; 


X  -  Propogata  alaaanta  - 

lagranga; 

i*i*l:  X  Updata  aatriz  indaz 


and  X  End  of  tha  intagration 

toe 

Xn-8iza(P,l); 

anargy-(.5*(P(347,5)+P(373,5))*8UB(P(348;372,5)))*ta  X  Trapazoidal  int. 
PlOO-anargy/tf  X  Povar  araragad  otar  antira  orbit 

pbar-anargy/3240  X  Povar  avaraga  ovar  thruat  pariod 


dv-(v-1.5*pi)*180/pi 
dl-(l-pi/2)*180/pi 
dRAAB- (RAAN-RAANo) * 180/pi 
da«(a-42162862)/1000 


X  Changa  in  arguziant  of  parigaa 
X  Changa  in  I 
X  Changa  in  RAAH 
X  Changa  in  a 
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82501 


X  Chang*  in  a 


alaB(:,l)adat(:.l); 
alaaCl.l)— 180; 
dat(l.l)«*'i80; 

subplot (2, 2,1) 
plot(dat(: ,1) ,dat(: ,S)) 
ylabalCdalta  argunant  of  parigaa') 

subplot (2, 2, 2) 

plot(dat(: ,l).forcas(:,8),*7' ,dat(: ,1) ,forcas(: .7) , *r* ,dat(: .1) ,forc*s(: ,8) 
ylabal  ( *  Accalaration  (■/s*c'‘2) ' ) 
subplot (2. 2, 3) 

plot(*l*n(: .1) .currant (: ,3) ,*!*■(: ,l),curr*nt(: ,4) , *r* ,al*B(: .1) ,cuxT*nt(: , 
5).'g’) 

ylabal  ClL/n  (A-n/kg)') 
subplot (2, 2, 4) 

plot(dat(:,l),P(:,l),’:*,dat(;,l),P(:,2),*y*,dat(:,l),P(:,3),'r*,dat(:,l),P 

(;.4),*g*.dat<:,l),P(:,S).*b») 

ylabaK’Spacific  povar  (U/kg)*) 

Xdisp(’Hlt  any  kay  to  continua*) 

Xpausa 

subplot (2, 2,1) 
plot(aloa( : , 1) ,alan( : ,4)) 
ylabal ( ’Inclination’ ) 

subplot (2, 2, 2) 
plot(*l«i( : , 1) ,*laa( : ,5)) 
ylabal C’RAAM’) 

subplot (2, 2, 3) 

plot (alaaC : . 1 ) . alanC : ,2) ) 

ylabal (’saai-Mj or  axis’) 

subplot(2,2,4) 

plot  (alaaiC : ,  1)  ,*lasi( :  ,3)  ) 

ylabal ( ’ accantricity ’ ) 
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The  script  lagrangea.a  is  called  by  Bolniya.a.  It  evaluates  Lagrange's  planetary 
equations  and  integrates  the  orbital  elements  used  in  the  Molniya  stationkeeping  study. 

X  This  file  integrates  the  orbital  eleaents  according  to  Lagrange's 
X  planetary  equations 

adot*2*e*sin(T) / (n*sqrt  < l-e“2) ) *Fr/aass«2*a*sqrt ( l-e“2)/ (n*r) eFu/aass ; 
edot«sqrt ( l-e‘2) esinCv) / (n*a) *Fr/aASS«sqrt ( l-e‘2) / (n*a'‘2*e) • (a‘2* ( i-e‘2) /r 
-r) *Fu/aass ; 

Idot*  (  (r*cos  )  /  (n*aea*sqrt  ( 1  -e‘2)  )  )  *F3/Bass ; 

RAAMdot* ( (r*sin («-»v) ) / (n*a*a*sqrt ( l*e*2)*sin(I ) ) ) •FS/aass ; 
vdotl«->(sqrt(l>o‘2)/(n*a*e))*(cos(T)*Fr/aaas) ;  X  1st  tera 

«dot3*-r*(l/tan(I))*sin(v-»T)*F3/(n*a*a*sqrt(l-e“2)«aaas) ;  X  2nd  tera 

vdot2«(sqrt(l-e"2)/(n*a*e))*(l-»l/(l-*^e*cos(v)))*sin(T)*Fu/aass;  X  3rd  tera 

«dot««dotl'»«dot2-»Hdot3:  X  All  3  teras  of  eqtn  13 

dM*n- (2*r/a- ( l-e“2) *cos (v) /e) *¥t/ (aass*n*a) 

-(l-e“2)*(l+r/(a*(l-e'2)))*sin(v)*Fu/(n*a*e*aass) ;  X  Eqtn  17 

X  Euler  integration 
a>a<»adot*t8 ; 
e«e«edot*ts ; 

N-dM*ts<»N; 

«aa4ud0t*ts : 

I«I+Idot*ts; 

IUAN-IUAN<»RAANdot*ts ; 


X  Debug  data  collection 

debug(i, :)>Ci  t  v*180/pi  edotl  sdot2  «dot3  sdot] ; 
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Fuaction  k«pl*r.B  solves  Kepler’s  equation.  It  is  called  by  Bolniya.a. 


function  Eakeplerdl.Eo.e) 

£■£0-  (£o-e«sln(£o)  -N)  /  ( l**e*coa  (£0)  ) : 

X  This  function  it  utod  to  ittrativoly  find  £  froa  M  and  a. 

X  £0  ia  a  guaaa  of  £.  Rapaat  the  function  until  £0  convargaa  to  £. 


Function  rotate  .■  generates  a  rotation  matrix  given  two  rotation  angles.  It  is  called 
by  inc.a,  raan.a,  and  aolniya.a. 

function  [R]  ■  rot  (phi.laa) 

R«Ccoa(laa)«coa(phi)  -ain(laa)  ■*coa(laa)*ain(phi) 
ain(laa)*cot(phi)  coa(laa)  -ain(laa)*ain(phi) 
ain(phi)  0  coa(phi)]; 

X  Thia  function  calculatea  a  rotation  aatriz  baaed  on  the  two  rotation 
X  anglea  phi  and  laa 
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Script  d«cay*l.B  was  used  to  study  the  decaying  solution  to  the  rendezvous  study. 


X  This  is  the  file  used  to 
X  rendezvous  problea 
rhor«2.8e-8; 
rhoc«2700; 

ssicalO; 

r-6S78000; 

n>l.i834e-03; 

alf adot-7 . 2924620S6e-S ; 

zdo>0 ; ydo«-3*zo*n/2 : 

tf«.99*2*pi/n 


investigate  the  decaying  solution  of  the 

X  Resistivity  (Ohn-B)  [rho.R] 

X  Conductor  aass  density  (kg/B"3) 

X  Crho.c] 

X  Conductor  nass  ratio  Db/b_c] 

X  Radius  of  target  (■) 

X  Mean  notion  (rad/soc) 

X  Angular  velocity  of  the  Earth's 
X  rotation  (rad/sec) 

X  Calculate  initial  velocities  as  a 
X  function  of  initial  position 
X  tf  is  the  tine  period  in  which  the 
X  rendezvous  is  to  be  performed. 

X  See  the  file  for  the  polynomial 
X  solution  for  an  explanation  of  the 
X  .99  value. 


B«[0  0  -2.8282e-0S]: 


X  B  field  vector 


ao«(al‘2)/4; 

s>al/2: 

bz«zo: 


X  The  ao  coefficient 
X  This  is  Cp]  in  chapter  5 
X  More  contsants 


z«zdo-»s*zo: 

by-yo; 

cy«ydo*s*yo ; 
j-0; 

for  ia0:30:tf 
t-i; 

x»exp(-s*t)e(bx+cx*t) ; 
xdaexp(-s*t)*('s*bx-s*cx*t-»cx) ; 

y»exp(-set)e(by+cy*t) ; 
yd«exp(-s*t)*(-s*by-secy*t*cy) ; 

f x»-2*n*yd-al*xd- (ao+3en*2) ex ; 
f ya2*n*xd-ai*yd-ao*y ; 
X(j*l,:)»Cx  y  xd  yd  t] ; 

u(j+l,;)aCfx  fy]; 


X  Integrates  the  equations.  The  time 
X  step  is  30  sec. 

X  x  soln 

X  y  soln 

X  radial  acceleration 
X  in-track  acceleration 
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ILl— fy/B(3) ; 

IL2-fx/B(3); 

IL3-0; 

Ty* (n- alf adot ) •r : 

Pi«-ILl*Ty*B(3) ; 

Pz>4*rhor*rhoc*aBc*lLl ‘2 : 
Py*4*rhor*rhoc*«Be*IL2“2 ; 

■»MiP»Pi*Px+Py; 

IL(J*1.;)-CIL1  IL2  IL3]; 

P(j*lf!)"CPi  Pi*Px  Py  auaP]; 

•nd 

subplot (2, 2.1) 

plot(X(: ,2)/1000,X(; .1)/1000. • . *) 
axisCC'l.l  1.1  -1.1  1.1]) 
zlsbsl(*y  (ka)') 
ylsbsl(*z  (ka)*) 

subplot (2. 2, 3) 
plot<dst(:,l),X(:.3),»:»,dst(:,l),X(:,4),*-.»); 
ylsbslCv.z  k  tr.y  (a/soc)*) 

sz«szis ; 

szisCCszd)  5256  sz(3:4)]): 
zlsboK’Tiao  (ssc)'): 

subplot(2.2,2) 

plot(dat(:,l),u(:,l),»;»,dst(:,l),u(;,2),»-.') 
Xsz>szis:mzis(CO  180  sz(1.3)  •z(l,4)]) 
ylsbolCf.x  ft  f.y  (a/sse)') 
zlsboK'Tiao  (ssc)') 

szaazis : 

szisCCuCl)  5256  sz(3:4)]); 


X  Ck^ps.r] 

X  Ckspps.thsts] 

X  Cka4^ps.3] 

X  t.rsl 

X  povsr  duo  to  inducsd  xoltogs 
X  radial  i‘*2B  povsr 
X  ia-track  i*2B  povsr 
X  total  povsr 


subplot (2, 2, 4) 

plot(dat(:,l),P(;.2),';',dat(:,l),P(:,3),»-.',dat(;,l).P(:,4)); 
ylaboK'Spscific  Povsr  (V/kg)*) 
azaazis ; 

azis([az(l)  5256  ax(3:4)]): 
zlabsK'Tias  (ssc)'); 
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■■slz«(P,l) ; 

•n«rgy*(.S*(P(l. :)'»?(■. :))«suB(P(2:a>l. :)))*30  X  Trapezoidal  integration 
Pmz-mz(P) 

Pbaraenergy/tf  X  Average  power 


Script  poly.z  was  used  to  study  the  polynomial  solution  to  the  rendezvous  study. 

X  Thie  ie  the  file  used  to  study  the  polynoaial  solution  to  the  rendezvous 
X  and  docking  problea 


rhor«2.8e-8; 
rhoc»2700 ; 

r«6578000; 

n«l.i834e-03; 

alf adot-7 . 292462056e-5 ; 

zdo*0 ; ydo»-3*xo*n/2 ; 

tf».99*2*pi/n 


B-CO  0  -2.8282e-0S]; 


j-0; 


X  Resistivity  (Oha-a)  [rho.R] 

X  Conductor  aass  density  (kg/a*3) 

X  Crho.cl 

X  Conductor  aass  ratio  Diii/B.c] 

X  Radius  of  target  vehicle  (a) 

X  Mean  action  (rad/sec) 

X  Angular  velocity  of  Earth’s 
X  rotation  (rad/sec) 

X  Initial  conditions  based  on  zo 
X  and  yo 

X  This  is  the  tiae  span  you  want  to 
X  perfora  the  rendezvous  in.  99X  of 
X  the  orbital  period  is  used  so 
X  coaparison  can  be  aade  to  a  2  bum 
X  rendezvous,  like  the  ezai^le  in 
X  Kaplan  (pp  114*115).  If  the  exact 
X  orbital  period  is  used,  Kaplan’s 
X  eqtns  3.54  and  3.55  becoae 
X  singular. 

X  B  field  vector.  It  is  constant  for 
X  the  assumptions  aade  in  Ch  5. 


ax^xo;  X  Constants  for  the  x  soln 

bx«xdo ; 
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cx«-3*xdo/tf -6*xo/tf “2 ; 
dx«3*xdo/tf ‘2*8*xo/tf “3 ; 
•x*-xdo/ tf “3*  3*xo/ tf “ 4 ; 

»y*yo; 
byaydo ; 

cy»-3*ydo/tf -6*yo/tf “2 ; 
dy*3*ydo/tf “2*8*yo/tf ‘3 ; 
•y«-ydo/tf“3-3*yo/tf‘4; 

for  i«0:30:tf 


t-i; 

x«xx+bx*t+cx*t “2+dx*t  *3*«x*t ‘4 ; 
xd«bx-»2*ex*t+3*dx*t“2+4*«x*t“3 ; 
xdd*2*cX'»8*dx*t-»  12*«x«t  *2 ; 

y«xy*by*t+cy*t‘*2-»dy*t ‘‘34«y*t  “4 ; 
yd»by+2*cy*t+3*dy*t '‘2*4*«y*t‘‘3 ; 
ydd»2*cy*8*dy*t  ♦  12*oy*t  “2 ; 

X(j+l,:)*Di:  y  xd  yd  t]; 


fx«xdd-2*ii*yd-3*n“2*x ; 
fy»ydd*2*n*xd; 

u(j+l,:)»Cfx  fy]; 

ILl— fy/B(3) ; 

IL2-fx/B(3) ; 

IL3-0; 

vya (n-alf adot ) *r ; 

Pi— ILl*Ty*B(3) ; 
Pxa4*rhor*rhoc*BBc«ILl‘2 ; 
PyM*rhor*rhoc*BBc*IL2*2 ; 
luaPaPi-^Px-^Py; 


X  Constants  for  tbs  y  soln 


X  Intagratas  tha  aquations,  lota  tha 
X  tiaa  stap  is  30  sac. 


X  X  soln 


X  y  soln 


X  Collact  tha  position  and  valocity 
X  data 

X  Coi^uta  tha  accalarationa 


X  Ckappa.r] 

X  Ck^pa.thata] 

X  Ck^px-33 
X  T.ral 

X  poaar  dua  to  inducad  voltaga 
X  radial  i*2R  powar 
X  in- track  i*^  povar 
X  total  povar 


IL(j+l,:)«[ILl  IL2  IL3]; 
P(j+l.:)»CPi  Pi+Px  Py  su*P] ; 
dat(J+l,l)»t; 
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•nd 


subplot (2, 2,1) 

plot(X(:,2)/1000,X(:,l)/1000,*.‘) 
szis(C-l.l  1.1  -1.1  1.1]) 
zlsbsK'y  (tai)*) 
ylsbslCz  (In)') 

subplot(2.2,3) 

plot (date : , 1) .X( ; .3) , • : * ,dat ( : , 1) .X< : .4) , •- . • ) ; 

ylaboK'v.z  t  t.y  (a/ssc)') 

az«azis; 

azisCCazd)  S2S6  az(3:4)]): 
zlabalCTias  (sac)*); 

subplot (2, 2, 2) 

plot(dat(;,l),u(;,l),':*,dat(:,l),u(:,2),’-.*) 
Xaz«azis;azis(CO  180  az(l,3)  az(1.4)]) 
ylabalCf.z  ft  f.y  (a/sac)') 
zlabalCTiM  (sac)') 

az«azis ; 

azis(Caz(l)  5256  az(3:4)]); 


subplot (2, 2 ,4) 

plot(dat(;,l),P(:.2),*;',dat(:,l),P(:,3).'-.',dat(:,l),P(:,4)); 
ylabslCSpaciflc  Povar  (H/kg)') 
az*azis ; 

azis(Caz(l)  5256  az(3:4)]); 
zlabalCTina  (sac)'); 

■•siza(P,l) ; 

anargy«(.5*(P(l, :)'>P(B,:))<»suB(P(2:a-l,:)))a30  X  Tr^azoidal  intagratiou 
Pmz-mz(P) 

Pbar«anargy/tf 


X  Avaraga  powar 


Apptndix  B.  Survey  of  Conducting  Metals 


Equation  (47)  was  derived  in  section  2.6.3.  The  last  term  includes  the  product 
PkPo  where  is  the  resistivity  of  the  conducting  material  in  Ohm-meters  and  Pe  is  the 
conductor  mass  density  in  kilograms/meter*.  To  minimize  the  last  term  of  Equation  (47), 
it  is  desirable  to  use  a  conductor  with  a  low  pj,pc  product.  To  find  such  a  conductor  we 
compile  p«  and  Pe  values  for  common  conducting  metals  (3:1324-1325)  and  compute  PmPe- 


Metal 

Pm 

(Ohm-m  xlO"*) 

Pc 

(kg/m*  X  10*) 

PmPc 

(Ohm-kg/m*  X  10"“) 

Aluminum 

2.66 

2.7 

7.2 

Copper 

1.68 

8.9 

15.0 

Gold 

2.42 

19.3 

46.7 

Silver 

1.63 

10.5 

17.1 

Aluminum  is  the  worst  conductor  in  terms  of  resistivity,  but  its  low  mass  density  makes  it 
the  best  overall  conductor  for  electrodynamic  propulsion. 
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